NASA- TM-1 10608 


y / 



A NUMERICAL STUDY OF FUNDAMENTAL SHOCK 

NOISE MECHANISMS 4. y / 


■MMM. eirtfi 

IV.0S tUWlUTMW 


A Dissertation 

Presented to the Faculty of the Graduate School 
of Cornell University 

in Partial Fulfillment of the Requirements for the Degree of 

Doctor of Philosophy 


(NASA— TM— 110608) A NUMERICAL STUOY 
OF FUNDAMENTAL SHOCK NOISE 
MECHANISMS Ph.D. Thesis - Cornell 
Univ. (NASA. Langley Research 
Center) 193 p 

G3/71 0049641 


by 


N95-2 7908 
Unc 1 as 


Kristine R. Meadows 
May 1995 



© Kristine R. Meadows 1995 
ALL RIGHTS RESERVED 



A NUMERICAL STUDY OF FUNDAMENTAL SHOCK NOISE MECHANISMS 


Kristine R. Meadows, Ph.D. 
Cornell University 1995 


The results of this thesis demonstrate that direct numerical simulation can predict 
sound generation in unsteady aerodynamic flows containing shock waves. Shock 
waves can be significant sources of sound in high speed jet flows, on helicopter 
blades, and in supersonic combustion inlets. Direct computation of sound permits 
the prediction of noise levels in the preliminary design stage and can be used 
as a tool to focus experimental studies, thereby reducing cost and increasing the 
probability of a successfully quiet product in less time. 

Direct simulation of sound generation in shocked flows is challenging because of 
the disparity in amplitude between the acoustic waves and shocks. These challenges 
are met by the implementation of a high-order accurate Essentially Non- Oscillatory 
(ENO) scheme which maintains high accuracy in smooth regions of the flow to 
minimize numerical dissipation of the acoustic waves while maintaining sufficient 
numerical dissipation at the shock for stability. 

This thesis reveals and investigates two mechanisms fundamental to sound gen- 
eration by shocked flows: shock motion and shock deformation. Shock motion is 
modeled by the interaction of a sound wave with a shock. During the interaction, 



the shock wave begins to move and the sound pressure is amplified as the wave 
passes through the shock. The numerical approach presented in this thesis is vali- 
dated by the comparison of results obtained in a quasi-one dimensional simulation 
with hnear theory. Analysis of the perturbation energy demonstrated for the first 
time that acoustic energy is generated by the interaction. 

Shock deformation is investigated by the numerical simulation of a ring vor- 
tex interacting with a shock. This interaction models the passage of turbulent 
structures through the shock wave. The simulation demonstrates that both acous- 
tic waves and contact surfaces are generated downstream during the interaction. 
Analysis demonstrates that the acoustic wave spreads cylindrically, that the sound 
intensity is highly directional, and that the sound pressure level increases signif- 
icantly with increasing shock strength. The effect of shock strength on sound 
pressure level is consistent with experimental observations of shock noise, indicat- 
ing that the interaction of a ring vortex with a shock wave correctly models a 
dominant mechanism of shock noise generation. 
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Chapter 1 


Introduction 


The motivation behind this research is to establish a better understanding of the 
physical nature of sound generation in shocked flows. Knowledge gained by in- 
creased understanding of the primary mechanisms involved in shock noise genera- 
tion has application to a variety of real aerodynamic problems. In helicopter and 
tilt-rotor applications, shock waves may form on the blade surfaces, and vortices 
from a preceding blade may interact with the shock resulting in impulsive noise. 
In supersonic engine inlets, combustion instabilities may result in shock oscillation 
leading to what is commonly known as “buzz”. Even subsonic transport aircraft 
often operate at conditions where the flow over the wing is transonic, and un- 
steadiness in the flow over the wing can cause the shock to oscillate and generate 
noise. 

The primary application in mind during the course of this research is the sound 
generation in supersonic jets. When jet engines operate at supercritical nozzle 
pressure ratios, shocks may form in the jet plume and turbulence interacting with 
the shock waves generates h,gh amplitude, broad-band noise, typically known as 


1 
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“broadband shock noise” or “shock-associated noise”. The expected shock noise 
during the climb-to-cruise operating condition of a supersonic civil transport may 
ruin its chance of ever being put in production. Shock noise is an important 
design issue because of the effects on community noise, aircraft interior noise, and 
structural fatigue. 

Research towards the understanding of noise generation in jet flows has been 
ongoing for over 40 years. Morley [1] in an investigation of sound intensity in the 
far field of turbulent jets, showed that the sound power is proportional to about 
the eighth power of the jet velocity. Lighthill [2] in his pioneering work on jet noise 
theory, provided a theoretical basis for the eighth power law noticed by Morley, 
and provided the basis for an understanding of other jet noise phenomena such 
as convective amplification. Lilley [5] provides a nice review of classical jet noise 
theory and the related experiments. There are two recent review papers specific 
to noise generation in supersonic jets which summarize significant contributions 
to the understanding of jet noise. The first, written from an experimentalist’s 
perspective, is the review article by Seiner [3]. The most recent, written from a 
theoretician’s perspective, is the article by Tam [4], The complicated nature of 
the supersonic jet flow makes development of a comprehensive theory difficult, and 
therefore the theories currently available for the prediction of shock noise in jets 
are largely empirical in nature. 

The approach taken in this research is to compute directly the sound generated 
by shock waves in supersonic jets. Advances in computer hardware resulting in 
more computational speed, affordable memory, and parallel architectures combined 
with advances in software resulting in more efficient, accurate algorithms and better 
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networking have made this type of approach currently feasible for relatively simple, 
two dimensional problems. 

Direct computation of aerodynamicaJly generated sound makes sense because 
sound is inherently a component of a fluid flow field. Therefore, the basic equations 
governing sound are the same as those governing fluid flow. For Newtonian fluid 
flows, these equations are the Navier Stokes equations, or, when viscous effects can 
be neglected, the Euler equations. As can be seen by Figure 1.1, the acoustic por- 
tion of a fluid flow field consists of small perturbations on an inviscid, compressible 
flow. The generation of sound waves, however, often involves viscous, nonlinear 
effects. Thus, an advantage of performing a direct simulation of the fully nonlinear 
equations, as opposed to using traditional acoustic methods which require that in- 
formation along some surface is provided through another means, (e.g. experiment 

or a separate computation), is that both the sound generation and propagation are 
computed in the same analysis. 

If current trends in increasing computational capacity continue, it will even- 
tually be possible to simulate directly sound generation in complicated, three- 
dimensional, unsteady flows of practical interest, such as supersonic jet flow. How- 
ever, direct computation is currently not feasible for routine study of supersonic 
jet flow because the flow contains far too many scales to be resolved in a reason- 
able computation. Therefore, a combined modeling- direct computation approach 

ken in this research. First, the complicated flow field of a supersonic jet is 
broken down into simple “model” problems, chosen to isolate mechanisms which 
are likely to generate sound. The modeling of the supersonic jet flow is described 
in Chapter 2. Chapter 2 begins by a description of features characteristic of su- 
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Figure 1.1: Some well-known equations of fluid mechanics. Note that the acoustic 
wave equation is a special case of the Navier Stokes equations. 
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person* jet noise and then proceeds to break the complicated flow down into 
several simpler problems which can currently be solved using reasonable computer 
resources. Chapter 2 also presents an analytical model for shock noise genera- 
tion. In tins analysis, Lighthill's equation is used to derive the monopole, dipole, 
and quadrupole source terms associated with an moving planar shock. The source 
terms are presented and then a Green’s function method is used to compute the 
far held sound associated with each of the source terms. It is found that all of the 
source terms can potentially contribute significantly to the far field sound. 

Chapter 3 addresses numerical issues and methods involved in the computation 
of flows with moving shocks. Computation of sound generation in shocked flows 
is chaUenging because high accuracy is required to maintain the acoustic portion 
of the solution, while dissipation is required at the shock to maintain stability. 
Most of the calculations presented in this thesis use a finite-volume implementation 
of an Essentially Non-Oscillatory (ENO) scheme which provides these important 
algorithmic features. Because this algorithm is the basis for most of the calculations 

of this research, and is not currently described in text-books, a brief description is 
included in this chapter. 

Once candidate sound generating mechanisms are identified and modeled, and 
the numerical methods are in place, the flow for each model problem can be com- 
puted and the acoustics extracted from the calculation. The model problems which 
are addressed in this work are those associated with shock noise generation. Chap- 
ter 4 describes the results obtained by an investigation of the first model problem: 
sound wave-shock wave interaction. Computations of the interaction of a sound 
wave with a shock wave in a quasi-one dimensional nozzle are presented which 
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show significant amplification of the sound wave pressure amplitude as a result 
of the interaction. The increase in amplitude of the acoustic pressure is shown 
to be a function of the Mach number upstream of the shock. Comparisons with 
linear theory are made for the small-disturbance calculations which validate the 
code. Results are also provided for the higher-amplitude cases. In addition, an 

energy analysis is performed which shows that acoustic energy is generated during 
sound-shock interaction. 

Chapter 5 presents results obtained by the investigation of the interaction of a 
vortex ring with a shock wave. This sound generating mechanism is more compli- 
cated than the plane sound-shock interaction case because the shock bends dur- 
ing the interaction, generating alternating compression-rarefaction-compression- 
rarefaction regions along the acoustic wave front. Flow parameters downstream of 
the shock are observed and from these observations the conclusion is drawn that 
both acoustic waves and contact surfaces result from the interaction. Analysis of 
the results shows that the acoustic wave spreads cylindrically, that the sound inten- 
sity is highly directional, and that the sound pressure level increases significantly 
with increasing shock strength. The effect of shock strength on sound pressure 
level is consistent with experimental observations of shock noise, indicating that 
the interaction of a ring vortex with a shock wave correctly models the physics of 
shock noise generation. 

Chapter 6 summarizes the significant findings of this work. 
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Chapter 2 


Modeling of Sound Generating 
Mechanisms in a Supersonic Jet 


2.1 Introduction 


The inherently complicated nature of jet flown makes both theoretical analysis 

and direct numerical simulation impractical for realistic flow Mach and Reynolds 

numbers. The structure of the flow field of a supersonic jet is very complicated, 

consisting of regions of laminar flow, turbulent flow, and transitional flow. Within 

the jet there exist a myriad of structures of disparate scales, such as turbulent 

eddies and shock waves, which make analysis of the fluid dynamics practically 

impossible for general flows. Highly accurate numerical simulation of these flows 

is impractical because of the disparity of the scales required to be resolved in the 
computation. 

In order to predict accurately the sound generated by complex jet flows, the 
essential elements of the jet fluid dynamics must be resolved in the computation. 


8 



9 



Figure 2.1: Schlieren photograph of an underexpanded jet. 

Experimental .Indies of noise generation in snpersonic jets have provided insight 
into what these essential elements may be. A result from an experimental investi- 
gation of sound generation in a supersonic axisymmetric jet is included in Figure 
2.1. This photographic image, provided by Dr. J.M. Seiner of the Jet Noise Labo- 
ratory at NASA Langley Research Center [1], was produced by a horizontal spark 
schlieren of an underexpanded supersonic jet ( ratio of exit pressure to ambient 
pressure is 3.05, and the je, design Mach number is 1.8) and confirms the complex 
nature of the je, Sow. As the flow exits from the nozzle, it rapidly becomes tur- 
bulent. For the underexpanded case shown here, shock waves are also present in 

the flow, and a Mach disk is evident approximately one jet diameter downstream 
of the nozzle. 

For the acoustician, the wave field outside the jet plume is of particular interest. 



10 


It is obvious that several types of waves are present, indicating that there are several 
types of sound generating mechanisms responsible. 

Two types of wave fields are intense enough to be readily visualized on the 
schlieren photograph. The first appears to emanate from a region close to the noz- 
zle exit near marker ‘A’ on the figure and radiate at approximately ±30 degrees 
from the jet axis. These waves propagate to the far field as acoustic waves, and are 
thought to be generated by supersonically convecting eddies which occur when the 
jet flow is heated and in unheated flows when the plume Mach number approaches 
two [2], The eddies create a form of sound known as “eddy Mach wave radia- 
tion”, which generally dominates the noise spectrum in directions of its dominant 
directivity. The primary directivity can be determined by computing the angle 
complementary to the Mach angle. (The Mach angle is determined from the eddy 
convection velocity: fi = sin -1 jj). 

The second type of wave outside the jet plume appears to emanate from the 
terminal locations of the shock waves in the mixing layer. One such wave field is 
readily seen emanating from a region in the proximity of marker ‘B’ in the figure. 
This wave field appears to be more omni-directional than the eddy Mach wave 
radiation, because the wave fronts appear to spread spherically from the point of 
generation. This second type of acoustic wave is believed to be generated by the 
passage of turbulence through the shock waves and is referred to as “broadband 
shock noise or “shock noise”. The broad-band nature of this type of noise is 
illustrated in Figure 2.2. 

Figure 2.2 shows the spectrum of a typical underexpanded supersonic jet at 
150 degrees from the jet axis. The spectrum is characterized by three types of 
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St = fD/Uj 

Figure 2.2: Spectrum of a typical supersonic underexpanded jet. 

phenomena which are highlighted in the figure. The figure shows that these three 
components of noise tend to fall into distinct frequency bands. At low Strouhal 
number (a dimensionless frequency formed by the product of the frequency and 
jet diameter divided by the fully expanded jet velocity) jet noise dominates. Jet 
noise is generated by turbulent mixing within the jet. The peak in jet noise, at 
a Strouhal number of about 0.1, is believed to be generated by the large scale 
structures within the shear layer and is typically referred to as eddy Mach wave 
radiation. Both large and small scales within the jet generate mixing noise. The 
mixing of the fine scales produces the background level of jet noise. 

The peak of the spectrum is known as jet screech. Screech typically falls be- 
tween the jet noise and the broadband shock noise in the spectrum. The accepted 
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explanation for screech was originally proposed by Powell [3], who suggested that 
screech is generated when the turbulence interacting with shock waves in the flow 
develops a self-sustained aeroacoustic feed-back loop. The resonant phenomenon 
which generates a high amplitude sound at a particular frequency and its har- 
monics, is maintained by the shedding of a disturbance at the nozzle exit by the 
passing of sound past the nozzle lip. Although Powell’s explanation of the screech 
phenomenon is useful in explaining general features of the phenomenon, there are 
features of jet screech which are still not understood [2]. 

The high Strouhal portion of the spectrum of Figure 2.2 is called broadband 
shock noise. Broadband shock noise is present in all shocked jet flows, and is due to 
the interaction of convecting disturbances with the shock waves in the jet plume. 
The broadband nature of this noise is due to the many scales of turbulent eddies 

in the flow. Experiments (e.g., [2], [5] ) show that most of the broadband shock 
noise is directed slightly upstream. 

Although jet engines are designed to be shock-free at design operating condi- 
tions, jet engines are often operated at off-design conditions. Thus, screech and 
broadband shock noise can contribute significantly to the sound being generated. 

2.2 Models 

Because jet flows are too complicated to be practically analyzed or directly simu- 
lated, the approach taken in this research is to identify elements essential to sound 
generation in supersonic jets, and analyze models of these mechanisms for insight 
into the sound generation processes. The models presented here are probably not 
the only models necessary for an in-depth understanding of sound generated by 
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supersonic jets. However, it is assumed that understanding isolated mechanisms 
will provide insight for reduction of noise levels within a real jet flow. 

The basis for the models is the experimental evidence described in the first 
section of this chapter. Experimental data indicate that noise is generated by 
different types of mechanisms within the jet. To gain insight into the fundamental 
mechanisms associated with sound generation, the complex flow of a supersonic 
jet is modeled as illustrated in Figure 2.3. The first class of sound generating 
mechanisms in the supersonic jet arises from the interaction of flow disturbances 
with the shock waves. The elements which model the sound being generated by 
the interaction of flow disturbances with shock waves are illustrated on the left 
side of the figure: (1) the interaction of a vortex with a shock wave, and (2) the 
interaction of a plane sound wave with a shock. Because this research focuses on 
the noise generated by shocked flows, these two elements will be pursued in detail 
in Chapters 4 and 5. 

Additional models of sound generating mechanisms which are not directly re- 
lated to the presence of shock waves are illustrated on the right side of the figure. 
The model presented in the upper right of the figure is flow past a wavy wall. This 
models the component of sound generated by large scale structures within the jet 
plume. Clearly, one would expect Mach waves to be generated by flow past a wavy 
wall, and in the model, these Mach waves are analogous to the eddy Mach wave 
radiation from supersonic jets. The model presented in the lower right of Figure 2.3 
is presented to model the interaction of turbulent structures within the jet, which 
is believed to be responsible for the background noise (“jet noise”) of supersonic 
jet flows. 
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Figure 2.3: Fundamental sound generating mechanisms of supersonic jet flow. 

2.3 Sound Generation by Shock Motion: 

Evaluation of Lighthill’s Source Term and 
Far Field Sound 


This section presents an analysis of the sound generation by shock motion in the 
context of Lighthill’s acoustic analogy [6]. 

2.3.1 Background 

Lighthill combined the equations governing the conservation of mass and momen- 
tum into the form of an inhomogeneous wave equation. A brief derivation of 
Lighthill’s equation follows. 

The physical law governing the conservation of mass written in indicia! form is: 


dp d[puj) _ 

8t + dxj V 


3 = 1,2,3 


( 2 . 1 ) 
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where p is the fluid density, uj is the j — th component of velocity, t is time, x{ 
is the coordinate in the i — th direction, and Q is the mass per unit time per unit 
volume injected into the fluid. 

The conservation of momentum equations are written: 
dui duj dPu 

p it + + Fi i ' i = 1 ' 2 - 3 < 2 - 2 ) 

where Fi is an externally applied force, the stress tensor = p8ij~a X] where 8{j is 
the Kroneker delta, and c x ;y is the viscous stress tensor, <r x j — , where 

e 'i = an d F an d A are the dynamic viscosity and second coefficient 

of viscosity, respectively. 

Taking the inner product of u, with Eqn. 2.1 and adding it to Eqn. 2.2, one 
obtains: 


8 £ u i + ej^) + ^ i ^ Fi + niQ 

OXj UXj 


dt 


( 2 . 3 ) 


Differentiating Eqn. 2.1 with time, subtracting the divergence of Eqn. 2.3, and 
adding and subtracting one obtains LighthilPs equation: 


gV SQ ajFi + ujQ) 

dG 00 dt dxi dxidxj 

where the LighthiU stress tensor is defined: 


(24) 


Tij = puiUj - (Tij + pSij - clepSij (2.5) 

The terms on the right hand side of 2.4 are the acoustic source terms. These 
terms represent sound generated by unsteady mass addition (^2), unsteady forces 
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' dxl nonlinear viscous effects, turbulence, and nonlinear propagation ( puiuj ), 
and non-isentropic effects such as shock waves and heat addition ((p — pc^Sij). 
Effects of viscosity on the sound generation are represented by a{j. 

2.3.2 Source Terms 


Here, Lighthill’s equation will be used to analyze the sound generated by shock 
oscillation. To simplify the analysis, dissipation effects are neglected, and only 
the velocity component normal to the shock, u, is considered. Since there are no 
applied forces or unsteady mass addition, the Lighthill source term simplifies to 


T n = puu + p - pc^ 


( 2 . 6 ) 


The source term of the Lighthill equation requires that the second partial deriva- 
tive of T n be taken with respect to x. It is beneficial to use generalized functions 
because ordinary derivatives do not exist across the shock. The definition of a 
generalized derivative in one dimension is [8]: 


Bfjx) 

dx 


9f{x ) 

dx 


+ A fS(x - z,) 


(2.7) 


where £ represents the generalized derivative operator, ^ is the ordinary deriva- 
tive, and Af represents the jump in the function / at the discontinuity located at 
* = x t . 

Applying 2.7 twice to obtain , 


dT 


li 


dx 


+ AT n S(x - x t ) 


8 2 T] 

dx 2 


n 


d^T ii d dT 1 

+ fc( ATll6 ( X ~ **)) + A (-^ i K( a: ~ *.) (2.8) 
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Figure 2.4: Schematic of shock oscillation used in the determination of Lighthill’s 
source term. 

Thus, the Lighthill source term for sound generated by a moving shock has com- 
ponents which resemble shearing stresses (quadrupole), unsteady forces (dipole), 
and mass addition (monopole), as represented by the first, second and third terms 
on the right and side of Equation 2.8, respectively. 

In classical acoustics, where sound generation is considered in an ambient 
medium, the sound pressure associated with monopoles, dipoles and quadrupoles 
is proportional to the second, third, and fourth power of the source Mach number, 
respectively. Therefore, the monopole source dominates at low Mach numbers, 
and the quadrupole source dominates at high Mach numbers. However, for sound 
generation at transonic speeds, each component may contribute significantly to the 
sound generation. It is a purpose of this section to analytically determine which 
term dominates the shock noise generation. 

To determine the behavior of the monopole, dipole and quadrupole terms as 
functions of Mach number, consider a planar shock, set in motion by an upstream 
disturbance in Mach number, as illustrated in Figure 2.4. The Mach number 
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upstream of the shock is 


M\ - Mo — M a cos(u >t — kx ) (2.9) 

where M 0 is the undisturbed upstream Mach number, M„ is the maximum ampli- 
tude of the disturbance Mach number, u is the frequency of the disturbance, and 
k is its wavenumber. 

Disturbances in pressure, density, and sound speed upstream of the shock are 
neglected because these quantities are proportional to the square of the upstream 
disturbance Mach number. 

Monopole Source Term 


The monopole term on the right hand side of Equation 2.8 is 




Considering for now only A(^ i ), 




<9Ti 


li 


dx 


dT n 


dx 


( 2 . 10 ) 




c\M. 


d 


2 + P2 - P2COO] - — [pic\ Ml +Pi- pic\ 


(4 Mi - cl)^ + 2 P2 c 2 M 2 


dx 
d( c 2 M 2 ) 

dx 


2 1 

ooJ 


, dp 2 0 2 X, 9 Ml 


dx 

( 2 . 11 ) 


since pi,pij and c\ are, to first order, constant. 

To evaluate Equation 2.11 in terms of the upstream density, sound speed and 
Mach number, the well-known shock jump relations (see, e.g. [7]) are used. These 
relations are presented below for completeness. The second equality refers to the 
case of an ideal gas with the ratio of specific heats, ■y = 1.4. 
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P 2 _ 2 jMl - (7 - 1) _ 7Mf - 1 
Pi (7 + 1) 6 


( 2 . 12 ) 


P2 _ (7 + 1 )M\ _ 6 Ml 

Pi ~ (7 - 1 )Ml + 2 - M[+ 5 


«2 _ (7 ~ 1 )M[ ± 2 _ M[+_5 
ui (7 + l)Mf 6 M\ 


To compute the Lighthill monopole source term, the derivatives of the down- 


stream pressure, density, and velocity are required. Taking the first derivative with 
respect to x, for the case of 7 = 1.4, one obtains: 


dp 2 

dx 



dMi 

dx 


5 2 dM 1 

3 P1C ' M '-dT 


5(c 2 M 2 ) ci k{Ml - 5) 6Mi 


dx 


Ml 


dx 


(2.15) 


(2.16) 


dp2 _ fin Mi dMi 

dx ~ Pl (5 + Ml) 2 dx 


(2.17) 


Q If 

where = —M a ksm(wt — kx). 

Substituting in Equations 2.12, 2.13, and 2.14 into 2.11, and working through 
the algebra one obtains: 


A( 


dT n 
dx ’ 


- 60 pic^ 


Mi dM\ 
(Mi + 5) 2 dx 
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Figure 2.5: Normalized monopole component of the Lighthill source term sound 
generated by a sinusoidally oscillating shock. 

„„ o MiM.k 

= ~ 60plC °° lMf^ 5)2 sm M ~ kx ) (2.18) 

Substituting M\ — Mo — M 9 cos(a;t — kx), and keeping only terms of order M S} 
one obtains: 


^ o MoM.k . , , x 

60pl C °° (Mq + 5)2 sm ^ ut ~ kx ^ ( 2 - ] 19 ) 

Equation 2.19, normalized by GOpic^k sin(u;t — A:x), is plotted as a function of 
Mo (for constant wavenumber) in Figure 2.5. 

Two cases are shown. The upper curve in the figure represents the source term 
when the ratio of perturbation Mach number to upstream Mach number, K is hdd 
fixed. The lower curve represents the source term for a fixed perturbation Mach 
number, M s . Clearly, the Lighthill monopole source term reaches a maximum 
value. To determine the exact location at which it maximizes, the first derivative 
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of Equation 2.18 with respect to Mq is determined and set equal to zero. The 
value of Mo for which this derivative is zero, holding M fixed, is found to be 
M) = \/| ~ 1-29. The value of M\ for which is zero, holding ^ 

fixed, is found to be Mq = y/b « 2.24. 

Thus, the monopole source term reaches a maximum when, for constant shock 
velocity amplitude, and constant disturbance wavenumber, Mo = y|. When 
the ratio of the shock velocity amplitude to the pre-shock Mach number is held 
constant, the source term maximizes at Mo = y/b for a fixed wavenumber. 

For an acoustic disturbance, the wavenumber varies with upstream Mach num- 
ber, 


k 


U) 

ci(l + Mi) 

U! 

ci(l + M 0 ) 


+ 


COs(fcx — Ujt) 

(1 + M 0 ) 2 


M s + O(Mg ) 2 


( 2 . 20 ) 


The wavenumber may be approximated by k = an d substituted into 

Equation 2.19 while maintaining a truncation error of order (M,) 2 . Furthermore, 
if the sound speed approaching infinity is taken to be the stagnation sound speed, 
the relationship between c\ and c Q 0 is: 


ci 




Mo cos(o;£ — kx) 
(5 + M 2 ) 



Thus, the monopole source term for k = Aj(Mq) is 


( 2 . 21 ) 


—GOpiCoouMoMg sin (o»t — kx) 

71 (M 2 + 5) 3 /2(i + Mo) 


( 2 . 22 ) 
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Figure 2.6: Normalized strength of monopole component of the Lighthill source 
term for a sinusoidally oscillating shock. Wavenumber k is a function of upstream 
Mach number and sound speed. 

Equation 2.22 normalized by -SOpi^uM, sin(u;< - kx ) is plotted in Figure 

2 . 6 . 

Dipole Term 

Consider now the dipole term of Equation 2.8. This term can be written: 

Q 

[A7ii£(® - x a )] = AT n 6'(x - x 3 ) (2.23) 

since the jump across the shock is not a function of x. Thus, the jump in is of 
interest. 

A Iii = P2u\ + P2 - P2C 2 od - [piul +P 1 - picle] (2.24) 

But, from the unsteady jump relations (see Appendix A for derivation), 

P2{u 2-u s ) = pi(ui-Ug) 
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P1U 2 {U2 - U s ) + p 2 ~ P\U\{ui - U s ) + 


P 1 


(2.25) 


one can solve for the change in pressure across the shock: 


Vi - Pi = Pi{ui - u a )(ui - u 2 ) 


(2.26) 


Substituting 2.26 into 2.24, and neglecting higher order terms in u„, one obtains: 


ATn ft: ~{p 2 - pi)c^, (2.27) 

Therefore, upon substitution of 2.13 into 2.27, and simplification, the following 
expression for the dipole source term is obtained: 


ATu = ~ 5p ' cl ° Osfrl) < 2 - 28 > 

Substituting M x = M 0 - M, cos (w< - kx ) into Equation 2.28 and neglecting 
terms of order M] and higher, one obtains: 


AT n = 


-5pic 


.2 

'OO 


(MS + 5) 


(MS - 1 ) 


12M,Mq cos(u >t — kx) 


(2.29) 


Mq + 5 

It is interesting that unlike the monopole term, the dipole term contains both 
zero-th and first order terms in M a . Another difference between the monopole 
and dipole terms is that the wavenumber affects both the amplitude and phase of 
the monopole term, while it affects only the phase of the dipole term. It is also 
interesting that the monopole term and the first order component of the dipole term 
are identical except for the phase and the multiplicative factor of the wavenumber. 

As illustrated in Figure 2.7, there are no local extrema for the dipole source 
term when the flow is supersonic. This can be verified analytically by taking the 
derivative of Equation 2.28 with respect to Mq as for the monopole case. 
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Figure 2.7: Normalized strength of dipole component of the Light hill source term 
for a sinusoidally oscillating shock. 

Quadrupole Term 

Finally, consider the quadrupole component of Eqn. 2.8. The quadrupole term 
requires that the second ordinary derivative of the Lighthill stress tensor component 
Til with respect to x be evaluated. Thus, 


d 2 T] 
8x 2 


li 


d du 9 dp 

d 2 u 


dp 


d P i 


+ _ _ c z UL 


= 2 ^ + iu: a* ax 


dx dx °°dx i 

dxdx +2p ^dx^ ~ C °°)fc2 + di 2 


(2.30) 


The evaluation of Eqn. 2.30 upstream of the shock is simple, since the deriva- 
tives in p , p, and c are high-order, and hence, neglected. The evaluation of Eqn. 
2.30 upstream is: 


d 2 T n 

dx 2 j 


= 2 p! c\ 



+ Mi 


d s M t 

dx s 


(2.31) 
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Substituting M\ — M 0 M 3 cos(u >t — kx ) and neglecting higher order terms 
M s , this expression reduces to 


m 


d 2 T \ 


11 


dx 2 


_ 10 p c 2 t.2 MoMgCosjut - kx) 

°° (5 + M 0 2 ) 


(2.32) 


For an acoustic upstream disturbance for which the wavenumber is a function 
of Mach number, this becomes: 

9 2 T\\ n 5 MqM 3 cos(u )t — kx) 

-wr p,w < 2 - 33 > 

The evaluation of the quadrupole term downstream of the shock is more compli- 
cated since the second derivatives with respect to a: of p 2 , u 2 , and p 2 are required. 
The second derivatives of the variables in Eqns. 2.12, 2.13, and 2.14 
d 2 p 2 60/3i 


are: 


dx 2 


(5 + Mj 2 y \ 

d 2 p 2 5 




dx 2 ~ 3 Cl Pl 


f dM,Y d s M t 

+ M ‘~air 


dMA 
dx ) 


d 2 u 2 


6 Ml p { dx ) ~ 


( 5 Mi - MY) 


3 *d* Mt 


(2.34) 


(2.35) 


(2.36) 


dx 2 6 M? \ \ dx I > dx s 

where ^ = -kM t sin(u >t - kx), and - k 2 M t cos(u it - kx). 

Substituting Eqns. 2.35, 2.34, and 2.36 into Eqn. 2.30 and algebraic manipu 
lation provides: 


d 2 T n 
dx 2 

= 2 pt cf 

2 

7 dMA 2 

Wr) + 


+ 12pi 4 

M t 4 +20 Mf - 


(5 + Mi 2 ) 3 


m 
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+ 


12/Ci c 


2 

OO 


- 5 Mj_ 

(5 + MS) 2 


d^Mi 

d x 2 


(2.37) 


Substituting Mi - M 0 - M, cos(u >t - kx ) and neglecting higher-order terms, 
one obtains: 


10pl c °° k2 (5 + ° M 2)2 ( M o - l)cos(u it - kx) (2.38) 

It is interesting to note that the quadruple source term contains only first order 
terms (and higher) in M, and goes to zero as M 0 -> I. 

Substituting in Eqn. 2.20 for k = k(M 0 ), the quadrupole term becomes: 


Mo M s 




(2.39) 


It is interesting to note that the quadrupole source term goes to zero as M 0 — > 1, 
and that the quadrupole and dipole source terms have the same phase. 


2»4 Evaluation of Acoustic Pressure 


The previous section provides expressions for the monopole, dipole and quadrupole 
source terms of Lighthill’s equation for shock motion in one dimension. In order to 
make statements regarding the relative importance of these components to the far 
field sound, a Green’s function approach is used to solve for the far field acoustic 
pressure. The primary assumptions in this analysis are that only the first order 
terms in perturbation Mach number M a are necessary to determine the relative 
importance of the monopole, dipole and quadrupole terms, that the observer is 
located in the far field, and that the shock is a finite disk in a plug flow, as 
illustrated in Figure 2.8. The flow outside the plug is ambient and is characterized 
by the density Poo and the sound speed c^. For the analysis presented here, 
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Figure 2.8: Schematic of shock disk in plug flow. Flow variables are functions 
of position along the axial coordinate. The shock is located at the center of the 
volume. 

contributions of the interface between the plug flow and ambient medium to the 
far field sound are neglected. 

2.4.1 Monopole 

In order to evaluate the far field sound resultant from the monopole source, a 
cylindrical coordinate system is set up with the origin at the shock disk center, 
as illustrated in Figure 2.9. The position of the source point along the disk is 
denoted by y and the position to the observer is denoted by x. The magnitude 
of the vector connecting the source and observer points R = x - y is denoted R, 
which is determined by geometrical arguments to be: 

R = sfr o + <r 2 — 2ro<r sin 0 cos(<^ — x) (2. 40) 

2 

Making the far field approximation — > 0 and employing the binomial expan- 

r o ^ 
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Figure 2.9: Coordinate system used in evaluation of surface integrals, 
sion, this simplifies to: 

R « ro[l - — sin#cos(<£ — x)] (2-41) 

TO 

where, as illustrated in Figure 2.9, ro is the distance from the observer to the center 
of the shock, a is the radial coordinate to the source point, 8 is the angle from the 
jet axis to the observer, 4> is the observer angle about the axis of symmetry, and 
X is the angle to the source point in the plane of the shock disk. Applying the 
Green’s function approach to solve for the far field acoustic pressure resulting from 
the monopole component of the Lighthill source term, the source term is multiplied 
by the free space Green’s function and integrated over the source volume V and 


time r: 
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n(S,t) = -± J j 60 W s^r-k yi )S( yi -j,,.) ‘ ± - ‘ ± V^l jVir 

(2.42) 

where y\ a — 0 is the shock position, and dV = dyi<rd<rdx is the incremental volume. 
Making the far field assumption such that ^ and integrating over y u 

™ {i ’ t) = iirro J 1 6 ° P1C ^(^+ 'If 5i "(“' r “ - < - 


A, r 


*/ c oo 


(2.43) 


where A. 3 is the shock area. Integrating over source time t, 

PM{£,t) = 47rro ~ R/coo — kyu)dA B 


4i rr 0 7 " ri 00 
-60pi<4 MoM„k 


4irr 0 (M$ + 5) 2 
R. 2 ir 

f f ■ / . / & sin 0 cos(^ — v) 

J J sin (ut - (r„ ( 2 . 44 ) 


0 0 


Integrating with the aid of expressions 3.719-2 and 3.715-13 in [9], the analytical 
expression for the monopole source pressure is: 


Pm(x, t) = 


-SOpic^kMoM, R, 


J\ (R a sm6 — ) sin(u>f - (2.45) 

V Cqo / Cno 


(M 0 2 + 5) 2 r 0 a; sin# 
where Ji is the first order Bessel function of the first kind. 

2.4.2 Dipole 

Again using the Green’s function approach to solve for the pressure due to the 
dipole component of the Lighthill source term: 
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, GOpic^MgMo , 

+ (m 0 2 + 5) 2 cos ( a;7 ' ~~ k yi) 6 (yi - yu) 


S(t — t + — ) 

c -^-dVdr 


R 


(2.46) 


Making the far field approximation £ « integrating the first term over r, and 
using the fact that [8]: 


J f(yi)f>'(yi - yu)dy\ = 


df 


dy i 


(2.47) 


vis 


the integral of the first term goes to zero. Making the far field approximation and 
integrating the second term over r: 


/-. .s If 60 pic 2 M,Mo , . 

Pd(x, t ) _ ^ J v _____ cos(a/(< - R/ Coo ) - kyi)S {y x - y u )dV (2. 


r o Jv (Mq + 5) 

Integrating over y\ using Eqn. 2.47: 

1 60pic*M,Mo t d 


48) 


PD(x,t ) = 


(M 0 2 + 5) 2 J dyi 
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/ ~d^ ~ ~ **0] l»i=o dA a 


= 1 GOpi^MsMo f 

r 0 (M 2 + 5) 2 J R l c °°)- k yi) 
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w ^ ere fjn = -cos ^> an d dA, = <rd<rdx. Integrating over the shock area, the 
expression for the dipole source contribution to the far field sound is: 

„3 
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It is interesting to note that 
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2.4.3 Quadrupole 


Using the Green’s function approach to determine the contribution of the down- 
stream part of the quadrupole source term to the far field sound, 


?«(*,<)=/ J 10 /» l 4* 2 - (M ; + ° 5)2 ,(ut - k Vl )dVd 


(2.52) 

where iZ = ^ fr g + y? + cr 2 — 2royi cos 6 — 2ro<r sin 0 cos x- Making the far field 

2 2 

approximation that ~ 0 and ~ 0, and using the binomial expansion, 
R ~ ro[l — ^ cos 0 — ^ sin 0 cos x]- Note that in integrating the quadrupole term, 
the integral is over the entire source volume, since there are no delta functions 
which limit the integral to the surface. This provides a challenge for this model 
because integrating over an infinite volume does not correspond to physical jet 
flows. Therefore, for the purposes of this analysis, the volume term is limited to 
a length —4 R s < yi < 4 R, along the axial coordinate. This range is selected 
to approximately correspond to average shock spacing in supersonic nozzles [10]. 
Therefore, for integration of the downstream quadrupole term, the integration 
limits are 0 < y\ < AR 3 . Attempts at obtaining an analytical expression for this 
integral were unsuccessful. Therefore, numerical integration is used to evaluate the 
quadrupole term. 

To simplify the comparisons and isolate the sound generated downstream on 
the sound field, only the region downstream of the shock is considered in the 
volume integration. For the purposes of determining the relative importance of 
the monopole, dipole, and quadrupole terms it is also assumed that the source is 
compact, i.e., the wavelength of the sound is large compared to the source size. The 
use of the compact assumption allows an analytical expression for the quadrupole 
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term to be obtained, since R in the argument of the cosine term in Equation 2.52 
reduces to r*o. 

Employing the compact source assumption, the analytical expressions for the 
monopole, dipole, and quadrupole terms are found to be: 




(Mq + 5 ) 2 r 0 ir 


PD(x,t) = p M (x,t) 
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(2.55) 


Figure 2.10 shows the normalized acoustic root-mean-square pressure of the 
monopole, dipole, and quadrupole terms for a particular case with an observer 
located 100 shock disk radii (ro = 100) away from the mean shock position at an 
angle of 0 — j from the jet axis. For ease of comparison, the root-mean-square 
pressure is normalized by the constant multiplicative factor of the monopole, dipole, 
and quadrupole pressure terms, namely: 

15 P1 4M,M 0 M, 

7rro(M^ + 5) 2 ' ' ' 

and the flow variables are normalized with respect to the sound speed c 0 c , the 
shock disk radius R a , and the upstream density p\. For these comparisons, the 
wavenumber, k, varies with the mean flow Mach number as prescribed in Eqn. 2.20. 
There are several features worth noting on this figure. First, the quadrupole term 
increases with mean flow Mach number more rapidly with increasing frequency. 
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Figure 2.10: Root-mean- square far field pressure of the monopole, dipole, and 
quadrupole terms in the Lighthill analysis of shock noise. Results for three fre- 
quencies are shown. 
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The dipole terms increase with mean flow Mach number. The variation in the 
monopole terms is an order of magnitude smaller than the variation in the dipole 
terms. Although difficult to see on this scale, the monopole terms peak and then 
decrease with increasing Mach number. The monopole terms dominate for the 
low frequencies over the entire Mach number range, and for the high frequency 
case at low Mach For the case w = 1, the quadrupole term begins to dominate 
at about Mo = 2.8. At low frequencies (for which the source is highly compact), 
the contribution of the quadrupole terms is negligible over the Mach number range 
considered in the analysis. The wavelength of the sound is 


27rc 00 \/5(l + Mo) 


(2.57) 


wy5 + Mq 

Thus, low frequencies correspond to highly compact sources. For example, for 
u — 1, the wavelength over the range of 1 < Mo < 3 is 11.5 < A < 15. The range of 
wavelengths are not much larger than the source size (4 R 3 = 4), and the compact 
assumption begins to break down. On the other hand, the range of wavelength for 
a source frequency of u> = 0.1 is 115 < A < 150, and the compact assumption is 
valid. 


Figure 2.11 shows the effect of observer position on the root-mean-square pres- 
sure of the monopole, dipole, and quadrupole terms. The results shown are for an 
observer located at 100 shock disk radii away from the mean shock position, at an- 


gles of 6 = 0, f , and f from the jet axis. The monopole and quadrupole terms are 
unaffected by the change in observer position; Pm at all observer positions shown 
is seen to change only modestly about a value of approximately 0.016, and pq is 
seen to increase with mean flow Mach number from 0 at Mo — 1 to approximately 
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Mean Flow Mach Number, M 0 


Figure 2.11: Root-mean-square far field pressure of the monopole, dipole, and 
quadrupole terms in the LighthiU analysis of shock noise. Results for three observer 
angles are shown. 
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0.003 at Mo = 3.8. The dipole component is affected by the observer position. The 
dipole term is stronger at observer angles closely aligned to the jet axis, which is 
consistent with directivity patterns typical of dipoles. In addition, the dipole term 
reduces to the monopole term at an observer position of 6 — 7r/2. 

2.5 Conclusions 

Detailed analysis and direct numerical simulation of the sound generated by realis- 
tic supersonic jet flows are currently impractical. One way to gain insight into the 
sound generating mechanisms of such a complicated flow structure is to model the 
elements considered to be essential to the sound-generation process. Four essential 
elements were presented: Flow past a wavy wall, vortex-vortex interaction, sound- 
shock interaction and vortex-shock interaction. Because the focus of this work is 
on noise generation in shocked flows, the last two of these elements flows will be 
studied numerically in Chapters 4 and 5. 

This chapter presented an analysis of the Lighthill source terms associated with 
a model of an oscillating shock. The results presented for the acoustic pressure 
are accurate to first order in the perturbation Mach number, and are valid only 
for an observer in the far field. The figures presented are valid for situations 
in which the source is compact. The analysis shows that the monopole, dipole, 
and quadrupole terms are all potentially important in shock noise generation. The 
term which dominates in a particular situation depends upon the observer position, 
frequency, and mean flow Mach number. For the cases tested here, the monopole 
term dominates for low frequencies over the entire range of Mach numbers studied. 
Both the monopole and dipole terms are significantly larger than the quadrupole 
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term at low frequency. At high frequencies, the quadrupole term dominates at high 
Mach number. 
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Chapter 3 


Numerical Methods and Issues 


3.1 Introduction 

This chapter addresses numerical issues concerning the computation of sound gen- 
eration in shocked flows. The direct simulation of sound generation in shocked 
flows is challenging because high accuracy is required to resolve the acoustic por- 
tion of the solution, while dissipation is required to maintain stability at the shock. 
Because shock generated sound places high demands on the algorithm, care must 
be exercised in the selection of the scheme used in the simulations. Two algorithms 
are used in the results presented in this thesis. The first is the popular MacCor- 
mack scheme which is second order accurate in time and in space. It is used as 
a baseline in scheme comparisons because of its computational efficiency and sim- 
plicity. MacCormack’s scheme is relatively inexpensive, but Gibb’s oscillations can 
occur in the vicinity of the moving shock, even with added artificial dissipation. 
MacCormack s scheme, as it is implemented in this work, is described in Section 
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3.2.1. The Essentially Non-Oscillatory scheme is used in most of the calculations 
presented in this thesis. This class of schemes is relatively new, and is not yet 
described in text books. It was chosen because it allows for high-order accuracy 
in smooth regions of the flow while minimizing oscillations around discontinuities 
in the solution by the use of adaptive stenciling. Second, third, fourth, and fifth 
order accurate implementations of the ENO scheme were used during the course 
of this research. The scheme is described in Section 3.2.2. 

Section 3.3 deals with numerical issues regarding the calculation of shocked 
flows. The Gibbs’ phenomenon produced by attempting to resolve a discontinuity 
on a finite mesh is well known, and is not described in this thesis. There is an 
additional oscillatory phenomenon which manifests itself when the shock wave is 
moving slowly relative to the mesh. This additional phenomenon has only recently 
been observed and described [2], and has only very recently been quantified for 
high-order schemes [3]. Much of the content of this section has been published by 
the author in [3], but is included here for completeness, with permission from the 
publisher. 

It is well known that a high degree of accuracy is required to perform acoustic 
calculations. High accuracy can be obtained by using many cells (or grid points) 
in a low-order accurate scheme, or by using fewer cells (or grid points) in a higher- 
order accurate scheme. A method of determining which approach is most eco- 
nomical is outlined in Section 3.4. The result shows that the most economical 
approach depends on the degree of accuracy required. For acoustic calculations 
where the degree of accuracy required is on the order of 10” 7 , the 3rd order ENO 
scheme proves to be most economical. Hence, this is the scheme used in most of 
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the calculations of Chapters 4 and 5.. 

Finally, Section 3.5 presents information regarding the extraction of acoustic 
phenomena from aerodynamic flows. The separation of acoustics from aerodynam- 
ics cannot be performed in the near field, and thus the properties of acoustic waves 
in the far field are given. The implications of these properties on setting up and 
running a calculation for sound generation are then described. Some methods of 
analyzing the acoustics are provided. 

3.2 Algorithms 

3.2.1 MacCormack’s Scheme 

The classical MacCormack’s scheme is employed in this research because it is ef- 
ficient and has been used extensively in computational aeroacoustics. (e.g. [17], 
[19], [18]). MacCormack’s scheme is described in many textbooks (e.g. [21], [22]), 
but for completeness the algorithm is briefly outlined here. 

Consider the system of hyperbolic equations: 


dU dF n 


(3.1) 


where U is a vector of J components and F is a J component flux function of U. A 
simple example of such a system of equations is the Euler equations. MacCormack’s 
scheme approximates the solution to Eqn. 3.1 by a two step predictor-corrector 
technique. For each component of the U vector, the predictor step is: 


or 1 =vt- - p n 


(3.2) 
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and the corrector step is: 

ur^wj+ur-^Ff-FU) (3.3) 

where the superscript » represents the n“ time level, the subscript j represents 
the j‘ k ceU, and F represents the flux function F evaluated at the predictor value 

u. 

MacCormack’s scheme is second order accurate in space and in time. Analysis 
of the stability of MacCormack’s scheme, with discussion about the amplification 
and dispersion factors can be found in textbooks such as [21] and [22], 

MacCormack’s scheme as presented in Eqns. 3.2 and 3.3 experiences difficulty 
maintaining stability when applied to shocked flows. Therefore, additional artificial 
dissipation is required. In this algorithm, artificial viscosity of the Jameson type 

[20] is added to the right hand side of the predictor step. The corrected numerical 
flux is defined by 


where 


*3+1 -Ff-F 


AV 


(3.4) 


FAV '•+i/2( £, '+l - u <) - 4?i/2( p i+2 - 3£/j+i + 3 Vi - £/,-_,) (3.5) 

and 

4+1/2 = K^ 2 \\u\ + c)imax[vi,vi+i\ 

4+1/2 = m ®*[0, /f (4) (|u| -f- c) t - — 4+1/2] 

where 

v _ I Pi+\ — 2 Pi -f- Pi_i| 

P«+ 1 + 2 pi + 


(3.7) 
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where p is the pressure, u is the flow velocity, and c is the sound speed. In the 
calculations presented here, the coefficients are defined K = I anc j /f(4) _ _l_ 

4 256 

The pressure term in is generally of second order in smooth regions of 
flow, the e( 4 ) term dominates, and the artificial dissipation is fourth order. Near 
discontinuities, however, the pressure term reduces to zeroth order, the e ( 2 ) term 
dominates, and the artificial dissipation is second order. 

3.2.2 ENO Scheme 

The class of ENO schemes is relatively new [7], and is chosen because high accuracy 
is achieved in smooth portions of the flow, while spurious oscillations around flow 
discontinuities such as shocks are bounded. A brief discussion of the properties of 
ENO schemes is given below, but interested readers are referred to the references 
[7] and [8] for details. 

An ENO solution operator E h is r th order accurate in the sense of local trun- 
cation error 

E h U n = U n+1 + 0{h T+1 ) (3.8) 

where U(x ) is the sufficiently smooth exact solution, and h is the cell spacing. 
The distinguishing property of ENO schemes is that spurious oscillations near 
discontinuities in U are bounded. For the one-dimensional scalar case, this can be 
written: 

TV{E k U ) = TV{U) + 0(h 1+9 ) (3.9) 

for some q > 0, where TV(U ) is the total variation of U as a function of x, defined 
as TV(U) = 53 ,■ \Ui+! — Ui\ for a discrete solution to a scalar conservation law. 
Bounding oscillations near discontinuities is accomplished in ENO schemes 
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through the use of an adaptive stencil. Because the stencil used for discretiz- 
ing the differential equations adapts to the solution, schemes based on the ENO 
property may be thought of as adaptive filters or non-linear algorithms. 

ENO schemes minimize numerical oscillations around discontinuities by using 
data from the smoothest part of the flow. At each cell, a searching algorithm 
determines which portion of the surrounding flow is smoothest. The stencil span- 
ning this portion of the flow is then used to construct a higher order accurate, 
conservative interpolation to determine the variables at the cell interfaces. 

In this particular finite-volume implementation, the interpolation operator is 
applied to the cell averaged characteristic variables, and the accuracy in space and 
in time is third order. Time integration is accomplished by the use of a third order, 
three stage Runga-Kutta scheme discussed by Shu [8]. The algorithm is applied 
to the conservation law form of the equations so that shocks are captured in the 
computations and no shock fitting is required to enforce the Rankine-Hugoniot 
jump relations across shocks which appear in the solution. 

3.2.3 Stencil Biasing Parameters 

ENO schemes can achieve high-order accuracy in smooth regions and capture 
shocks without oscillations by the use of adaptive stenciling. As the ENO schemes 
were originally presented in [7], the stencil shifts freely due to any detection of a 
numerical gradient. The direction of the stencil shift is determined by the mag- 
nitudes of the neighboring divided differences. The stencil shifts away from the 
larger differences. However, a loss of accuracy can occur when this freely adap- 
tive algorithm is used [10] [9], Shu [11] has suggested that the stencil be biased 
towards a preferred stencil, the one that makes the scheme stable in the sense of 
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linear stability analysis. The stencil is allowed to shift only when one neighboring 
difference is larger than the other by some factor. This factor will be referred to 
as the “bias” parameter. 

This biasing can be accomplished by implementing a factor in the search stencil. 
To explain how biasing is used to affect the stencil, let Af denote the operator that 
yields the k-th forward difference on a stencil of k + 1 cells with a left-most index 
i, which is defined recursively by: 

Aju = ui-f-i - Hi 

A k { u = A^u - A^u, A: = 2,3,... (3.10) 

The algorithm begins by setting j^(i) = i. This one-cell stencil results in a piecewise 
constant reconstruction that is spatially first order accurate. In order to choose 
j k+1 (i), k = 1 , ...r — 1 , the two stencils considered as candidates are those obtained 
by annexing a cell to the left or right of the previously determined cell. The stencil 
which is selected is the one in which the k — th difference is smaller in magnitude: 

f O - * 1 ’ <a R \A k R u\ 

3 + (*) = < 

{ otherwise. 

where A and A j^u are the k-th differences obtained by annexing the cell to the 
left or right of the previously determined stencil, respectively, (ctlj&r) = (1>^) or 
(o-, 1) for biasing to the left or right, respectively, with & > 1. 

Even when a bias parameter is used, there may be a loss of accuracy when all 
the numerical gradients in a region are small, but some are orders of magnitude 
larger than others. Atkins [12] has suggested the use of another parameter, which 
serves as a threshold, to force the shift to the preferred stencil whenever neighboring 
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differences are small, regardless of their relative magnitudes. This parameter will 
be referred to as the ^threshold” parameter. 

The threshold parameter is implemented as 

if|A£u|<e and|Ajju| < e then j k+1 (i) =j£ +1 (i) ( 3 .H) 

where e is a small parameter and j|(i) identifies the stencil obtained by annexing 
the k ~th cell in the linearly stable direction. 

Research regarding the appropriate stencil biasing parameters is on-going [16]. 

However, the results presented in this thesis used the definitions cited in the articles 
of [12] and [11], 

3.2.4 ENO Flux Computation 

ENO schemes of Godonov type rely upon the solution of the Riemann problem to 
calculate numerical fluxes. Two methods of computing the fluxes across the cell 

interfaces are evaluated in this thesis. The first is due to Roe [13], and the second 
is due to Osher [14]. 

Consider a system of hyperbolic equations, 

dU dF 

«r + fe =0 (312) 

where U is a vector of J components, and F is a J component differentiable function 
of U. 

Roe’s approximate Riemann solver determines the change in flux by finding a 
mean Jacobian matrix A which satisfies: 


A F = AAU 


(3.13) 
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where A represents the difference between any two states in solution space. The 
matrix A is required to have a complete set of right eigenvectors and reduce to 
the exact Jacobian when the states to the left and the right are equal: A(U, U) = 
9u(U). If Ay), r^\ and Swj are the jth eigenvalue, right eigenvector and inner 
product of the jth left eigenvector with A U, respectively, the Eqn. 3.13 can be 
written: 

AF = £ f^X (j )Swj (3.14) 

y= i 

The flux at the cell interface is written: 


f - 2 { Fl + Fr) - | Ayj | SwjrM (3.15) 

3 

where F L and F R are the values of the fluxes to the left and right of the interface, 
respectively. Roe gives the expressions for A, f and S Wj for the Euler equations in 

[13]- 

Osher’s approximate Riemann solver computes fluxes in state space rather than 
physical space. The flux difference between the left and right states is written: 


A F = 


L 


Ur dF 
u L 8U dU 


(3.16) 


where the integral is evaluated along an arbitrary path T in state space. 
The flux at the cell interface is given by: 


f = + f m |A w |r % (3.17) 

The evaluation of the integral in Eqn. 3.17 requires knowledge of the states along 
each subpath T(j) and any sonic states that occur. Osher and Solomon [14] solve 
for these states explicitly for the Euler equations. 
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3.3 Numerical Error Generated by a Slowly 
Moving Shock in a Duct 

3.3.1 Introduction 

The numerical treatment of unsteady shocks is challenging. In addition to the usual 
concerns of stability and accuracy, there are conflicting requirements regarding the 
calculation of the shock region. On one hand, it is desirable to resolve the shock 
crisply by minimizing the smearing effect of artificial dissipation at the shock. On 
the other hand, some artificial dissipation is required at the shock to minimize or 
eliminate the oscillations which occur when attempting to resolve a discontinuity 
on a finite mesh. The Gibbs phenomenon produced by attempting to resolve a 
discontinuity on a finite mesh is well known, has been discussed widely in relation 
to steady shock calculations, and will not be addressed here. There is an additional 
oscillatory phenomenon which can manifest itself in the computation of slowly 
moving shocks. This additional oscillation resulting from the unsteady nature of 
slowly moving shocks has only recently been discussed in the literature [1] [2] [4], 
and is investigated here to determine the nature of the spurious oscillations and 
the effect that these spurious numerical oscillations have on computing sound. 

In order to simplify the analysis and better isolate difficulties in the numerical 
calculations, only one dimensional and quasi-one dimensional flows will be treated 
here. Thus, vorticity waves will not appear and emphasis will be on predicting 
the acoustic and entropy waves. The model problem to be investigated is a shock 
moving at a constant velocity in a one-dimensional flow field. 

Spurious oscillations in unsteady computations of slowly moving shocks have 
been described by Woodward and Collela [1], who observed the oscillations in com- 
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putations of high-pressure ratio (p 2 / Pl > 10 5 ) shocks. These osciUations appear 
when the speed of the shock, relative to the mesh, is small compared to the max- 
imum flow speed at the shock. Woodward and Collela suggested that additional 
numerical dissipation be added to the scheme at the shock, and explained that 
the reason for this spurious numerical behavior is that the shock transition layer 
alternates between being thick and thin as it passes through the mesh. Roberts [2] 
has noted that the oscillation phenomenon is not observed in discontinuous solu- 
tions of scalar equations. Roberts explains the oscillation in terms of the discrete 
shock structure and shows that among first order flux difference splitting schemes, 
Osher’s approximate Riemann solver provides the smallest oscillations because the 
unsteady nature of the numerical shock structure in state space most closely ap- 
proximates the true shock structure. Roberts observed the spurious oscillation 
phenomenon in calculations for shock pressure ratios as low as 1.2 in calculations 
using first order flux-difference splitting schemes. Lindquist and Giles [4] observed 
spurious oscillations in computations using a Jameson style Runge-Kutta scheme 
with blended second- and fourth-difference artificial dissipation and also with a van- 
Leer flux vector splitting algorithm for shock pressure ratios 1.5 < P 2 /P 1 <2.1 
They described the oscillations in terms of the changing shock shape as the shock 
traverses the computational mesh, and found, as did Woodward and Collela, that 
the spurious oscillations could be reduced by smearing the shock over more com- 
putational cells. 

Spurious oscillations have been observed during the course of this research in 
computations using the classical MacCormack scheme, a recently developed high- 
order accurate Essentially Non- Oscillatory (ENO) scheme, and a recent imple- 
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mentation of Jameson’s Runge-Kutta scheme which employs a symmetric Total 
Variation Diminishing (TVD) matrix dissipation. Because the interest here is to 
compute sound generated by shocks and shock-fluid interactions, artificial dissipa- 
tion is not explicitly increased over the solution domain because of the deleterious 
effect on the generation and propagation of sound in the solution. This has created 
some difficulties because the spurious oscillations are preserved in the high-order 
accurate flow computations. 

ENO schemes of Godunov type rely upon the solution of the Riemann problem 
for the calculation of numerical fluxes. The effect that two Riemann solvers have on 
the spurious oscillations will be described in this thesis. In addition, because ENO 
schemes use an adaptive stencil to reduce spurious oscillations, modifications of this 
adaptive procedure will be examined. Finally, the effects of shock pressure ratio, 
Courant number, grid spacing, and shock speed on the amplitude and frequency 
of the spurious numerical oscillations will be described. 

3.3.2 Analysis: Exact Solution 

The governing equations for the inviscid, compressible flow in a constant area duct 
are assumed to be the one-dimensional Euler equations, 

dU 8F 

-m + d^ = 0 ( 318 ) 

where U is the vector of conserved variables [p, pu, pe] T , and F is the flux vector, 
[pu, pu 2 4- p, (pe + p)u} T . Standard notation is used; p is density, u is velocity, 
e is total energy per unit mass, and p is pressure. 

Consider these equations along the duct length 0 < x < L for t > 0, with the 
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initial condition: 


U(x, 0) = 


Ui x < 


(3.19) 


where the constant states 1 and 2 represent the flow upstream and downstream of 
a shock, respectively, and x a is the shock position. 

If these equations are integrated along the duct, one obtains: 
d f L 

I I J rt .ti 1 EV TT { T ±\\ Tl/TT/n . \ \ - 

(3.20) 


UJx + F(U(L,t))-F(U( 0 ,t)) = 0 


Oae solution is the trivial case, U\ = U 2 . The non-trivial solution to this 
equation with the initial condition 3.19 is a flow Held with a shock moving at 
constant velocity, u a , which satisfies 


-^(^ 2 ) F(U\ ) — u a {U2 — Ui) (3.21) 

Eqn. 3.21 follows directly from Eqn. 3.20 for constant area states on either side of 
the shock. The solution U\ = U 2 also satisfies Eqn. 3.21. 

3.3.3 Results 


The unsteady, compressible, inviscid flow in the duct is solved by numerical mte- 
gration of the one dimensional Euler equations. All variables in the supersonic flow 
field at the inflow boundary are prescribed. The static pressure at the downstream 
(subsonic) boundary is prescribed. Variables are normalised by the duct length, 
stagnation pressure and stagnation sound speed. Control over the shock velocity 
is obtained by making a transformation u = u - u„ where u. is the prescribed 

shock velocity, so that a positive shock velocity moves the shock to the left of the 
computational domain. 

Computations have been performed over a range of shock pressure ratios and 
shock speeds, but in the interest of brevity, only one typical calculation is shown 
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here. The calculations were performed on a 512-cell grid at a Courant number of 
1. Unless noted otherwise, calculations are performed using Roe’s flux solver and 
both the bias parameter and threshold parameter are ‘on’, meaning that the stencil 
is biased towards the preferred one, and a threshold limit is set. The values of the 
biasing parameters used in these calculations, referring to Equations 3.2.3 and 3.11, 
are a — 2 and e = 10 -3 . Figure 3.1 illustrates the spurious oscillations observed 
with the 3rd order ENO scheme for a case in which the pressure ratio across the 
shock is 10.333 and the shock is moving to the left at a speed of 0.05. Figure 
3.1 shows the pressure and entropy distributions in the duct after the shock has 
moved 15 percent of the duct length. Entropy is measured by the quantity s = . 

Although there were no oscillations at the shock in the initial shock position, 
once the shock begins to move, spurious waves develop in the flow solution. The 
oscillation is seen primarily in the entropy wave; the pressure wave is relatively 
unaffected. 

The spurious error is due to the discrete motion of the shock moving through 
the mesh. When the shock is located at a cell interface, it is extremely thin. As it 
moves through the cell interior, it smears out, weakens in strength, and entropy and 
pressure waves convect downstream. If the “shock passing frequency” is defined as 
the frequency associated with the shock passing through a cell, 


f shock 


U, 

Ax 


(3.22) 


where Ax is the grid spacing, and u, is the shock speed relative to the grid, 
the wavelengths associated with the pressure and entropy waves, A p and A,, are 
determined by 
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Figure 3.1: Pressure and Entropy as Functions of Distance along Duct. Shock 

Speed - .05, Shock Pressure Ratio = 10.33. Dashed line represents pressure. Solid 
line represents entropy. 
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where u 2 and a 2 are the downstream velocity, relative to the shock, and sound 
speed, respectively. When adequately resolved on the mesh, the spurious pressure 
and entropy waves measured in the numerical computations compare well with the 
wavelengths A p and A, described by the above equations. For example, the wave- 
length A„ of the entropy wave in Figure 3.1 is computed: A„ « .550 * .00195/0.05 « 
.021. Inspection of Figure 3.1 verifies this result. 

A senes of calculations, for the same p 2 /pi = 10.333 shock strength, illustrates 
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the effect of shock speed on the behavior of the spurious oscillations. These calcu- 
lations were performed on a 512-cell grid at a Courant number of 1.0. Figure 3.2 
shows the entropy distributions in the tube after the shock has moved to the left 
for a normalized time of 3.0 for shock speeds of 0.02, 0.05 and 0.15. For clarity, the 
entropy distributions are offset by a constant value of 0.1. The entropy upstream 
of the shock is 1.0. The wavelengths of the perturbations downstream of the shock 
are consistent with Equation 3.24. The long wavelength disturbances at the lowest 
shock speeds are only slightly damped downstream of the shock, while the short 
wavelength disturbances at high shock speeds are damped very quickly by the dis- 
sipation in the scheme. This, of course, explains why these disturbances are not 
seen when the shock speed through the mesh is comparable to the flow speed. 

Figure 3.3 shows the pressure distributions in the tube after the shock has 
moved to the left for a normalized time of 3.0 for shock speeds of 0.02, 0.05 and 0.15. 
For clarity, the pressure distributions are offset by a constant value of 0.05. The 
pressure upstream of the shock is approximately 0.027. Small perturbations are 
visible in the pressure distribution downstream of the shock, and the wavelengths 
of these perturbations are consistent with Equation 3.23. The long wavelength 
disturbances at the lowest shock speeds are only slightly damped downstream of the 
shock, while the short wavelength disturbances at high shock speeds are damped 
very quickly by the dissipation in the scheme. 

Effect of Shock Speed 

Figure 3.4 summarizes the effect of the shock velocity relative to the grid on the 
maximum amplitude of the spurious entropy. The ratio of the magnitude of maxi- 
mum zero-to-peak entropy error amplitude to the jump in entropy across the shock 
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Figure 3.2: Entropy as a Function of Distance along Duct and Shock Speed. Shock 
Speeds = .02, .05, 0.15; Shock Pressure Ratio = 10.33. 
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Figure 3.3: Pressure as a Function of Distance along Duct and Shock Speed. Shock 
Speeds = .02, .05, 0.15; Shock Pressure Ratio = 10.33. 
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is plotted as a function of the ratio of the shock velocity (normalized by upstream 
stagnation sound speed) to upstream Mach number. These calculations were per- 
formed for a shock moving to the right at shock speeds from 0 to 5.0 on a 256 
cell mesh. The spurious entropy amplitude is machine zero when the shock is sta- 
tionary relative to the grid, increases to a maximum when .005 < jj* < 0.1, and 
decreases as the ratio of ^ increases further. 

Figure 3.5 summarizes the effect of the shock velocity relative to the grid on 
the maximum amplitude of the spurious pressure. The ratio of the magnitude of 
maximum zero-to-peak pressure error amplitude to the shock strength is plotted 
as a function of the ratio of the shock velocity (normalized by upstream stagnation 
sound speed) to upstream Mach number. These calculations were performed for a 
shock moving to the right at shock speeds from 0 to 5.0 on a 256 cell mesh. The 
results show that the maximum pressure amplitude is relatively independent of the 
shock speed, and that the maximum error in the downstream pressure relative to 
the shock strength is less than 0.15 percent for all of the cases studied. 

Figure 3.6 shows the zero-to-peak pressure error at a location approximately 
50 cells downstream of the shock. This figure shows that for faster moving shocks, 
the amplitude of the pressure error is more rapidly damped. This is consistent 
with the observations made regarding Figure 3.3, which relates the wavelength of 
the spurious pressure with the shock speed. 

Because the amplitude of the error introduced by a slowly moving shock is 
manifest primarily in entropy, further discussion will focus on this flow variable. 
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u/M 

Figure 3.4: Effect of Shock Velocity on Spurious Entropy. Three shock strengths 
are shown. 

Effect of Shock Strength 

Figure 3.4 also illustrates that the magnitude of the maximum spurious entropy 
amplitude is a function of shock pressure ratio. As the shock strength increases, 
the magnitude of entropy error, relative to the static entropy jump, decreases. In 
the weak shock case, M = 1.1, the spurious oscillations are close to 100% of the 
static jump in entropy over a range of shock speeds. This is because for the weak 
shock cases, the entropy jump across the shock is very small. 

The absolute levels of maximum entropy error are actually orders of magnitude 
higher for the higher Mach number flows. For example, at = 0.04, the absolute 
levels of maximum entropy error are Ss max = 0.00033, 0.0039, and .097 for Mach 
numbers M = 1.1, 1.3, and 3.0, respectively. 
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Figure 3.5: Effect of Shock Velocity on Spurious Pressure. Three shock strengths 
are shown. 
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Figure 3.6: Effect of Shock Velocity on Spurious Pressure. The error is measured 
approximately 50 cells downstream of the shock. Three shock strengths are shown. 
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Effect of Courant Number 

A sequence of calculations was performed with Courant numbers between 0.1 and 
1.0 (the stability limit) to determine the effect of Courant number on the spurious 
entropy. Neither the amplitude nor the wavelength of the spurious oscillations is 
found to be sensitive to the Courant number. This is consistent with the results 
reported by Roberts for flux-difference splitting schemes [2], 

Effect of Stencil Biasing Parameters 

In this section the effects of the stencil biasing parameters on the spurious entropy 
are investigated. To illustrate the effects of these parameters on the algorithm, 
Figure 3.7 shows a space-time diagram of the stencil used in the ENO algorithm 
for four combinations of stencil biasing parameters: (a) threshold on, bias on, (b) 
threshold off, bias on, (c) threshold off, bias off, and (d) threshold on, bias off. For 
clarity, the space time diagrams are limited to the region near the shock, and the 
total number of cells in the computations was reduced to 256. The shock is initially 
located at a: = .5 and moves with a velocity of .01 to the left. The upstream Mach 
number is 1.3. For cases (a) and (d), the threshold is on, and a centered stencil 
is used for the majority of the computational space. The centered stencil is the 
linearly stable stencil in this case. Cells on either side of the shock use downwind 
or upwind stencils as appropriate. For cases (b) and (c), the threshold parameter 
is off, and there is a great deal more stencil shifting in the smooth regions of the 
flow. It is interesting to note that the stencils in these two figures foUow entropy 
and acoustic wave paths downstream of the shock. The stencil traces for these 
wave paths have slopes corresponding to acoustic and convecting wave speeds. 
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Figure 3.7: Effect of Stencil Biasing Parameter and Threshold Parameter on Stem 
cil. White represents a downwind stencil; Black represents an upwind stencil; Gray 
represents a centered stencil. Case(a): Biasing on, Threshold on; Case(b): Bias- 
ing on, Threshold off; Case(c): Biasing off, Threshold off; Case(d) Biasing off, 
Threshold on. 256 Cells 
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Effect of Mesh Spacing 

The effect of the mesh spacing is illustrated by comparing Figures 3.8 and 3.9. The 
results presented in F lg ures 3.8 and 3.9 are for a shock (upstream Mach number 
is 1.3) moving to the left at a speed of u, = 0.01. The difference in the results 
occurs because of the difference in discretization. A 256 cell mesh is used to 
obtain the results of Figure 3.8, but a 1024 cell mesh is used to obtain the results 
of Figure 3.9. The entropy distributions in Figures 3.9 and 3.8 are offset by a 
constant of .0075 and the numerical values of entropy are removed from the vertical 
axis for clarity. The difference between tick marks is .005. Refining the mesh 
reduces the magnitude of entropy oscillation, and shows the effect of the biasing 
parameters more distinctly. Each combination of biasing parameters has a unique 
spurious entropy pattern. Case (a) is highly oscillatory, with multiple frequencies 
per wavelength of the oscillation. Case (b) has the smallest amplitude of entropy 
peaks. Case (c) has the fewest peaks per period of spurious oscillation, but the 
amplitude of the largest peak is high. Case (d) has large entropy peaks as well 
as multiple frequencies per oscillation. Although all of the results show significant 
entropy error, the results obtained by biasing the stencil and turning the threshold 
off (case (b)) provide the lowest amplitude of entropy error. 

Another effect of the mesh spacing is the reduction in the wavelength of the 
spurious entropy. In Figure 3.8, the wavelength of spurious entropy computed by 
Eqn. 3.24 is % .34. Refining the mesh in Figure 3.9 reduces this wavelength to 
A, « .086. The number of points per wavelength of the spurious entropy wave is 
the same for both computations, since [ s constant. 

Figure 3. 10 is included to show the effect of the Osher flux solver on the spurious 
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Figure 3.8: Entropy as a Function of Duct Distance for Various Combinations of 
Biasing and Threshold Parameters. Case(a): Biasing on, Threshold on; Case(b): 

Biasing on, Threshold off; Case(c): Biasing off, Threshold off; Case(d) Biasing off, 
Threshold on. 256 Cells. 



Figure 3.9: Entropy as a Function of Duct Distance for Various Combinations of 
Biasing and Threshold Parameters. Case(a): Biasing on, Threshold on; Case(b): 

Biasing on, Threshold off; Case(c): Biasing off, Threshold off; Case(d) Biasing off, 
Threshold on. 1024 Cells. 
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Figure 3.10: Entropy as a function of Distance Along Duct Length for Roe and 
Osher Flux Solvers. Biasing on, Threshold off. 

entropy. The result is shown only for case (b) because other cases provided similar 
results. The effect of the Osher solver on the spurious entropy is to remove one 
of the peaks (seen in the Roe solver results) in the entropy wave each oscillation 
wavelength. 

This section has illustrated numerical difficulties in computing slowly mov- 
ing shocks. Although changes in flux computation, algorithmic parameters, and 
Courant number were made in an effort to remove these spurious entropy waves, 
they persisted. 

Calculations of supersonic jet noise will often have shocks moving slowly relative 
to the grid. Although the results presented thus far have shown that spurious 
entropy exists in calculations with slowly moving shocks, it has also been noted 
that spurious pressure waves are very small in amplitude. 



64 


3.3.4 Summary 

The computation of slowly moving shock waves produces spurious, numerical en- 
tropy. The spurious entropy is a function of the algorithm used in the calculation, 
and, as seen by the modifications made to the ENO scheme, even slight changes 
in a basic algorithm can produce marked changes in the structure of the spurious 
entropy. It is interesting to note that this phenomenon has been observed by the 
authors in implementing the MacCormack scheme, ENO schemes, and a matrix 
dissipation Runge-Kutta scheme, and by others using flux vector splitting schemes 
[4], flux-difference splitting schemes with Gudonov, Roe and Osher flux solvers [2j, 
and PPM methods [1J. Spurious entropy normalised by the entropy jump across 
the shock decreases with increasing shock strength and increasing shock velocity, 
but is insensitive to Conrant number. The amplitude of spurious entropy pertur- 
bations is relatively unaffected by the type of flux solve, used, but the Osher solver 
reduces the number of peaks in the spurious entropy wave form. 

Because the amphtude of the spurious entropy wave is a function of the shock 
speed relative to the grid, the obvious method of eliminating the spurious wave is to 
move the computational grid with the shock during the calculation. Although this 
is not difficult in one-dimensional problems, it is considered unreasonable for the 
multi-dimensional problems of practical interest, and was not implemented during 
the course of this research. Another approach to reducing the spurious entropy 
is to increase the dissipation of the algorithm, as suggested by Woodward and 
Collela. This was not implemented, because the added dissipation would affect 
the acoustic waves as well as the entropy. Acoustic pressure waves seem to be less 
sensitive to numerical errors generated by slowly moving shocks. 
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3.4 Economics of Higher-Order Schemes 

The desired accuracy in a computational solution determines whether it is more 
economical to use a second order scheme with many grid points or a high-order 
scheme with fewer grid points. Fignre 3.11 illustrates how the cost of implementing 
these algorithms, as measured by CPU seconds per time step, varies as a function 
of the solution error. This figure is constructed in the Mowing manner. The 
quasi -one dimensional Euler equations are solved for the nozzle problem described 
m detail in Chapter 4. Grid refinement studies for the 2nd, 3rd, and 4th order 
ENO algorithms are performed on a Cray-2 computer. For each successively re- 
fined spatial discretization, the error in Mach number is computed for a steady 
isentropic flow in a converging-diverging nozzle. In addition, CPU time per time 
step is measured for the solution on each mesh. The relationships of error and 
CPU time as functions of grid spacing are then combined to construct the figure. 
The symbols on Figure 3.11 represent numerical results from the grid refinement 
studies. The fines are extrapolations of the numerical results over the error range 
of interest. Because the temporal and spatial accuracy properties are equivalent 
for these algorithms, the trends presented in Figure 3.11 will be similar whether 
the accuracy computations are performed on steady or unsteady flows. The error 
in Figure 3.11 is measured in the Li norm, defined by: 

E \M(x n ) - M{ Xn )\ 

n — 1 

N (3-25) 

where M(x) is the exact value, M( x ) is the numerical approximation, and N is 
the number of discrete values in the numerical solution. 

Figure 3.11 shows that when an error of order HT 1 is acceptable, the second 
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Figure 3.11: Computational Time per Time Step as as Function of i, Error for 
2nd, 3rd, and 4th Order ENO Schemes. 

order scheme is more efficient than the fourth order scheme by almost an order of 

magnitude. If, on the other hand, an error on the order of only 1(T« is acceptable, 

then the higher order schemes are an order of magnitude more efficient. It is 

interesting to note that continually increasing the order of accuracy in a scheme 

does not necessarily result in a significant reduction in CPU time, even at very low 

acceptable levels of error. Figure 3.11 shows that for the ENO scheme, most of 

the benefit is realized in moving from second to third order in the range of error 
norms considered. 

It should be noted that the results of Figure 3.11 are particular to the quasi- 
one dimensional ENO algorithms used in this study. Figure 3.11 will not apply to 
the implementation of all schemes, because the relative cost of obtaining high-order 
accuracy for different algorithms will vary. However, the procedure for determining 
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the relationship of computer cost as a function of acceptable error for different 
algorithms will be the same. 

3.5 Extracting Acoustics from Aerodynamic 
Flow 

3.5.1 Acoustics Defined 

In order to be able to extract acoustics from an aerodynamic flow, it is important 
to define the characteristics of the acoustic waves. For the purposes of this paper, 
acoustic waves have the following characteristics: 

(1) . Acoustic waves propagate to the far field at the sun. of the local velocity and 
sound speeds. 

(2) . Acoustic waves are isentropic. 

(3) . Acoustic waves are small amplitude. 

(4) . Acoustic waves maintain constant power. Thus, the decay in pressure of acous- 
tic waves corresponds to geometric spreading losses. This implies that acoustic 
pressure decays proportional to J in three dimensional flows and -fc in two dimen- 
sional flows; there is no decay in pressure amplitude in one- dimension. (Note that 
these decay properties are valid when wave steepening is negligible. The spreading 
characteristics of nonlinear waves are discussed in Whitham’s text [24].) 

These characteristics of acoustic waves can be readily measured in the far field 
where the amplitudes of the acoustic disturbances dominate those of the “hy- 
drodynamic” (non-propagating) disturbances; however, the application of these 
definitions near the source become more difficult. In fact, separating the acoustic 
disturbances from the hydrodynamic disturbances in the acoustic near field is not 
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considered to be possible because these quantities do not exist independently [23], 

Both hydrodynamic and acoustic disturbances are heard by a near field observer, 

but because the hydrodynamic disturbances decay more rapidly with distance from 

the source than the acoustic disturbances, only the acoustic disturbances impact 
the far field. 

3.5.2 Requirements for Acoustic Calculations 

The characteristics of acoustic disturbances determine the requirements of the di- 
rect numerical simulations used: 

(1). Highly-accurate. The small amplitude of the acoustic waves relative to the 
underlying mean flow and the large distances the acoustic waves travel require 
that the numerical scheme used is high-order accurate and/or uses a grid fine 
enough to sufficiently resolve the acoustic waves. The disparity of length and time 
scales typically of importance in acoustics makes accuracy a real issue, since small 
scales are often important and cannot afford to be filtered from the computational 
solution. In addition, boundary conditions must be accurately imposed. 

(2). Large computational domain. Because acoustic disturbances are defined as 
the propagating portion of the unsteady pressure field, calculation of acoustics 
requires a large computational domain. The size of the computational domain 
must be large enough that the pressure field can be measured at least one acoustic 
wavelength away from the source. Acoustic wavelengths are often much larger than 
the source(s) which generate the sound. 

(3). Long time solution. Long time calculations are necessary to compute at least 
one period of the far field sound if the acoustic signal is periodic. For non-periodic 
signals, seven to ten cycles of the lowest frequency are required for reasonably accu- 
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rate spectral estimates. (For a good description of accuracy in spectral estimates, 
see Hardin’s text [25].) 
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Chapter 4 


Interaction of Sound with a 
Shock Wave 

4.1 Introduction 

In this chapter the interaction of a sound wave with a shock is considered. This 
study is meaningful because it models the planar oscillation of a shock wave, and 
shock waves in real aerodynamic flows are inherently unsteady. 

In the first three sections, shock capturing formulations of the MacCormack and 
ENO schemes are applied to the governing equations of fluid dynamics for time 
dependent, shocked flow through a nozzle. The emphasis of these computations 
is to predict numerically the amplification of sound by a shock and to compare 
the solutions provided by the different algorithms. Validation of numerical results 
is always important, but is particularly important for these calculations because, 
although high-order schemes ensure high-order accuracy in smooth regions of the 
solution, they necessarily reduce to first order at the shock. Hence, accuracy is 
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lost in the vicinity of the most important region of the solution: the sound source. 
Since a linear theory exists for this problem, numerical results are validated for the 
case of small amplitude acoustic waves incident on the shock. Much of the content 
in the first three sections was originally published by the author [1], but is included 
here for completeness, with permission from the publisher. 

Section 4.4 includes an analysis of the disturbance energy associated with a 
sound wave passing through a shock. This analysis provides insight into the source 
of disturbance energy generated at the shock. 

4.2 Analysis 

4.2.1 Linear Theory 

Within the context of linear theory, only entropy and acoustic waves may exist in a 
quasi-one dimensional flow [2], In a fully linear flow, these waves are independent. 
However, when the flow field contains nonlinear features such as shock waves, the 
nonlinearity acts as a coupling mechanism between the linear waves. Thus, the 
presence of shock waves in a flow field makes it possible for a sound wave incident 
upon a shock to suddenly change its amplitude and generate an entropy wave. 

The physical explanation for this phenomenon is as follows. Entropy is gener- 
ated across the shock. For a steady flow, this change in entropy is independent of 
time. However, when a sound wave impinges upon the shock, the shock begins to 
oscillate, the change in entropy is no longer constant, and a periodic entropy wave 
is generated which convects downstream at the local flow velocity. In addition, the 
impinging sound wave is amplified. 

Linearized analyses of the interaction of small disturbances with shock waves 
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have been made independently by Blokhintsev [3], Burgers [4], Moore [5] and Powell 
[6]. Landau and Lifschitz [7] report that the ratio of transmitted to incident sound 
waves determined by linear theory is: 

Sp2 = (Mi + 1) + 2( 7 - - 1) - (Mi + 1)[( 7 - 1 ]M[ + 21 

Spi (M 2 + 1) 2( 7 - 1)M 2 2 (M! 2 - 1) - (M 2 + 1)[( 7 - l)Mj 2 + 2] ^ 

where Mi is the Mach number upstream of the shock (henceforth called the pre- 
shock Mach number) , M 2 is the Mach number downstream of the shock, and 7 is 
the ratio of specific heats of the fluid. Equation 4.1 as well as a similar expression 
for the ratio of static pressures, jj are plotted on Figure 4.1 for the ratio of specific 
heats, 7 = 1.4. Note that Equation 4.1 predicts an amplification of the acoustic 
signal as it propagates through the shock wave for all pre-shock Mach numbers. 
This is not surprising, since the mean flow pressure also increases across a shock, 
but it is interesting that the ratio of the perturbation pressures, is not the 
same as that for static pressures, 

4.2.2 Riemann Analysis 

The results of the previous section are based on linear theory in which disturbance 
amplitudes are assumed to be small. To determine the range for which the linear 
result is valid, one may compute the ratio of disturbance pressures across the shock 
by considering the iterative solution to the Riemann problem. The Riemann anal- 
ysis is performed by considering the space-time diagram illustrated in Figure 4.2. 
Initially a steady shock wave separates states 1 and 4. At some time At, a distur- 
bance is introduced upstream of the shock which moves the shock and produces 
acoustic and entropy waves downstream. Knowing the initial states 1 and 4, and 
the incident perturbation amplitude, and utilizing the facts that acoustic waves 
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Figure 4.1: Ratios of Static and Perturbation Pressures as Functions of Pre-Shock 
Mach Number. 

are isentropic and entropy waves introduce no pressure perturbation, the ratio of 
the perturbed states 2 and 1 may be found by numerical iteration. For an incident 
sound wave of pressure perturbation 6p\ = epi sin wt, the ratio ^ is of interest 
and is compared with the linear theory result in Figure 4.3 for several perturbation 
amplitudes and shows excellent agreement for perturbation amplitudes less than 
e = 10 1 . (Results for e < 10 -2 are visually indistinguishable from the Blokhintsev 
line.) Even for perturbation amplitudes on the order of e = 1.0, there is only a 10 
per cent difference between the solutions at M = 3. 

Combining the expression for ^ with the Rankine-Hugoniot shock jump re- 
lations, one may calculate similar expressions for the fluctuations in density and 
entropy downstream, as well as the shock velocity. These relations are plotted in 
Figures 4.4 and 4.5. Figure 4.4 shows the relationship between the ratios of pertur- 
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Figure 4.2: Diagram of the quasi-steady Riemann problem analysis for sound-shock 
interaction. Bold line represents shock. Dashed line represents entropy wave. Lines 
between states 3 and 4 represent acoustic wave. 



Upstream Mach Number, M 


Figure 4.3: Pressure Perturbation as Function of Upstream Mach Number. 
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Pre-Shock Mach Number, M 


Figure 4.4: Perturbation Ratios as a Function of Upstream Mach Number. 

bation pressure, entropy and density perturbations across the shock. These ratios 
were computed using a perturbation amplitude of e = 10 -3 . This figure shows that 
the pressure fluctuation is most significant over a wide range of Mach numbers. The 
entropy fluctuations downstream of the shock are very small; the quantity ^ is on 
the order of e, except near M = 1 where As — > 0. Figure 4.5 shows the difference 
between the upstream Mach number and the shock Mach number as a function of 
the upstream Mach number for several perturbation amplitudes. The linear theory 
results are presented in fines with no symbols, while the results obtained by using 
to obtained by the Riemann analysis are presented with dots. As expected, the 
results obtained for small perturbations (e < 10 -1 ) are indistinguishable from the 
results of the Riemann analysis, while there is a significant difference in the results 
obtained for the large perturbation (e = 1.0). The shock motion increases with 
upstream Mach number and perturbation amplitude. 
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Upstream Mach Number, M 


Figure 4.5: Shock Speed Number as a Function of Upstream Mach Number. Solid 
lines correspond to e = 1.0; Long dashed line corresponds to e = 10 -1 ; Dashed line 
corresponds to e = 10 -2 . Dots on lines represent results obtained from Riemann 
analysis. 
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Figure 4.6: Nozzle Geometry for Mach Slope =1. 

4.3 Model Problem 

4.3.1 Governing Equations 


The equations governing the unsteady qnasi-one dimensional flow in a nozule 


are: 


8U 6F 
dt 


dx~' v (4.2) 

where V is the vector of conserved variables |„, p is the flux yector 


[puA, ( pu 2 + p)A, ( pe + p)uA} T , where A is the 


area of the nozzle, and Q is 


source vector due to the area variation, [0, pg, 0]*. Standard notation is used; 
P is density, u is velocity, e is total energy per unit mass, and p is pressure. 

4.3.2 Nozzle Shape 


Consider a nossle which is converging-diverging, and designed for a linear Mach 
number distribution when flow is isentropic and fully expanded to produce super- 


(L-^L 
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Table 4.1: Mach slopes and corresponding range of pre-shock Mach number. 


Mach Slope 

Pre- Shock 

a 

Mach No. 

i 

1.4 to 1.8 

2 

2.0 to 2.6 


sonic velocity in the diverging section. The Mach number at the nozzle inlet is 
M = .8 and varies as a function of distance along the nozzle length: 

M(x) = ax + .8, 0 < x < 1 (4 3) 

where x is the distance along the nozzle, normalized by the nozzle length. Two 
Mach slopes are used to facilitate consideration of a range of practically significant 
pre-shock Mach numbers and keep the shock close to the center of the nozzle 
diffuser. The values of Mach slope, a, used in the calculations presented here and 
the corresponding pre-shock Mach numbers are presented in Table 4.1. 

A sketch of a nozzle with a Mach slope of 1 is shown in Figure 4.6. 

4.3.3 Governing Equations and Boundary Conditions 

The unsteady, compressible, inviscid flow in a nozzle of varying cross sectional area 
is governed by the quasi-one dimensional Euler equations 4.2. These equations 
are solved in conservation form so that shocks are automatically captured. Total 
pressure and entropy are prescribed upstream of the shock at the nozzle inlet, and 
static pressure is prescribed downstream of the shock at the nozzle exit plane. Flow 
quantities are normalized by the upstream stagnation conditions and the nozzle 

length, L. Finite wave conditions developed by Atkins and Casper [8] are employed 
at the boundaries. 
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The first step in numerically modeling the interaction of a sound wave with a 
shock is to determine an accurate steady flow solution throughout the nozzle. To 
obtain sufficiently accurate unsteady results, the steady state residual is converged 
to several orders of magnitude smaller than the smallest perturbation amplitude 
to be investigated. Steady flows with residuals only one or two magnitudes smaller 
than the perturbation amplitude may introduce spurious entropy at the inflow 
boundary in the unsteady calculation. 

To induce unsteady flow through the nozzle, the inflow boundary condition is 
perturbed sinusoidally and isentropically. The pressure perturbation is prescribed 


8p = e p sinwt (4-4) 

where e and w are the normalized amplitude and frequency of the perturbation, and 
p and t are the normalized pressure at the inflow boundary and time. The pressure 
perturbation introduces an acoustic wave at the nozzle inlet which propagates 
downstream at a velocity equal to the sum of the local velocity and sound speed. 
A range of perturbation amplitudes is selected so that linear as well as nonlinear 
waves can be investigated. The time dependent features of the flow are computed, 
and the effect of the shock on the sound wave is observed. 

4.3.4 Algorithms 

In order to admit discontinuous solutions, shock capturing formulations of the Mac- 
Cormack and ENO algorithms described in Chapter 2 are employed for the com- 
putations. Therefore, shock fitting methods which explicitly invoke the Rankine- 
Hugoniot jump relations across a shock are not required. 
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4.4 Results 

4.4.1 Unsteady Calculations 

The results presented in this section are for the case of perturbation amplitude, 
e = .01, and a frequency, u> = 60, which corresponds to an approximate wave 
number of 5.25. The calculations are performed on a grid of 128 cells, which 
corresponds to approximately twenty-four points per wavelength. The nozzle back 
pressure is prescribed so that a shock appears at x = .6 in the steady solution. For 
a nozzle with a Mach slope a = 1, this corresponds to a pre-shock Mach number 
of 1.4. 

Figure 4.7 shows the perturbation pressure, density, and velocity in the nozzle 
at a normalized time of 7 r, and compares the results provided by the MacCormack 
and 2nd, 3rd, and 4th order ENO schemes. Perturbation quantities are determined 
by subtracting the mean flow quantities from their time dependent counterparts: 

6p(x,t) = p(x,t)-p(x, 0) 

6p(x t t) = p(x,t)-p(x, 0) (4.5) 

Su(x,t) = u(x,t) — u(a;,0) 

The following are observations of the numerical results: The MacCormack and 
3rd and 4th order ENO schemes do an excellent job of predicting the perturbation 
pressure upstream of the shock, where the flow is smooth. The second order ENO 
scheme, which behaves similar to a TVD scheme because it has only second order 
interpolation, has a slight leading phase error. At the shock, the ENO schemes do 
a better job of capturing the shock in fewer cells. The flow solutions differ more 
significantly downstream of the shock. The phase shift error is amplified and the 
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wave amplitude is dissipated in the 2nd order ENO results. The MacCormack 
solution is slightly damped relative to the 3rd and 4th order ENO solutions, and 
has spurious oscillations on portions of the perturbations close to the shock. The 
difference in the solution between the third and fourth order ENO schemes is not 
graphically perceptible on this scale. These results show that all the algorithms 
applied to this problem do a good job of computing the perturbation quantities 
in the nozzle. The 3rd and 4th order ENO schemes provide the best solutions 
because they have no oscillations at the shock and less damping downstream. The 
cost of running high-order ENO is higher than second order ENO or MacCormack, 
however, and the economics of high-order schemes is discussed in Section 2.4. 

4.4.2 Effect of Mach Number 

Sound amplification by a shock wave is highly dependent upon the pre-shock Mach 
number. To illustrate this numerically, Figures 4.8 and 4.9 show a sequence of 
pressure perturbations moving through the nozzle for pre-shock Mach numbers 
of 1.58 and 2.36. A small amplitude perturbation of e = 10~ 5 is introduced at 
the inflow boundary. The perturbation is seen to maintain its sinusoidal shape 
throughout the nozzle. The amplitude decreases as the flow expands in the nozzle 
upstream of the shock. At the shock, it is amplified as the flow is compressed, and 
continues to increase gradually as the flow is compressed further. It is clear that 
the amplification of the sound wave is much more significant at the higher Mach 


number. 
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Figure 4.8: Several Snapshots in Time of Pressure Perturbation Along Nozzle, 
e = 1(T 5 , Mx = 1.58. 



Figure 4.9: Several Snapshots in Time of Pressure Perturbation Along Nozzle, 
e = 1CT 5 , Ml = 2.36. 
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Table 4.2. Pre-shock Mach number and minimum cells per wavelength for the 
calculations presented in this paper. 


Pre-Shock 
Mach No. 

Minimum 
CeUs /Wavelength 

1.4 

26.92 

1.58 

26.18 

1.7 

25.75 

2.0 

23.52 

2.36 

22.50 

2.6 

22.13 


4.4.3 Comparison of Numerical Results with Linear 
Theory 

To quantify the effect of Mach number on the perturbation amplification and val- 
idate the computations, the numerically determined pressure perturbation ampli- 
tude ratio is compared in Figure 5 with the linear theory over a range of Mach 
numbers, and for two perturbation amplitudes: e = 10 -5 and e = 1CT 2 . The first 
amplitude is small enough to be well within the validity of the Unear theory. The 
second amphtude approaches the Umit of Unear theory vaUdity, particularly for the 
lower Mach numbers where the shock is weak and may undergo large excursions. 
The normahzed acoustic perturbation wavenumber at the inflow boundary is set 
to 4. Results for Mach numbers between 1.4 and 2.6 are presented here. The 
computations have from 22 to 27 ceUs per wavelength, depending on the particular 
cases being calculated. Table 4.2 below Usts the cells per wavelength for the cases 
presented in this paper. 

Figure 4.10 compares the numerical results with Unear theory. The 3rd order 
ENO scheme and Unear theory match extremely well for the very small amphtude 
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Pre-Shock Mach Number 

Figure 4.10: Pressure Perturbation Ratio as a Function of Pre-Shock Mach Num- 
ber. 

perturbation, e = 10 -5 . The differences between numerical solutions and linear 
theory become more pronounced at the higher perturbation amplitude, e = 10~ 2 . 
Some discrepancy between the linear theory and numerical results is not surpris- 
ing since nonlinearities may become important at this perturbation amplitude. 
Because the shock in the MacCormack solution is spread over several cells, deter- 
mination of the ratio ^ is more difficult, particularly when the shock is located 
close to the end of the nozzle. This is why no MacCormack result is shown for 

M - 2.6. 

4.4.4 Energy Analysis 

The previous sections have shown that according to linear theory and numerical 
computation, sound pressure is amplified as a sound wave passes through a shock. 
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This section will address the issue of whether 


or not acoustic energy is generated 


in the interaction process. 

To answer this question, consider Myers' energy coroUary [9J. Myers' energy 
corollary is an exact equation governrng the transport of energy associated with an 

arbitrary flow field. This coroUary aflows for a descrrption of the acoustic energy 
m general situations where the linear description of, 
with the highly nonlinear flow field associated with 
an approach is warranted here. 

For a one-dimensional flow in which 
lary reduces to: 


energy is inadequate. Clearly, 
sound-shock interaction, such 

viscous effects are negligible, Myers’ corol- 


oE dW 

dt ^ dx D (4.6) 

where the disturbance energy density is defined 

E 1 

A = 2^“ 2 - “») + Pou o ■ (no - ») - pT o{ , - „) + p(h _ K) (< ?) 

Where 4 is the cross sectional area of the duct, T is the absolute temperature, h is 
he enthalpy, a is entropy, p ls density, p is pressure, u is velocity. The subscript „ 
represents the undisturbed state. The disturbance energy flux is 
W 

A PU P ' lU,J >i h ~ h «~ T 0 (s - + mUo (T - r 0 )(a - ( 4 8 ) 

and the disturbance energy source is 


D _ 

A - (pu - pqu q )( s - s 0 ) • VT 0 - (s - sq)pouq • V(T - T 0 ) (4.9) 


For the one-dimensional 
modes are present (vorticity 


sound-shock interaction, only acoustic and entropy 
requires three dimensions). Thus, the energy density 
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can 


be divided into two types of energy, acoustic and entropy, as: 


1 \ dh 
Ea = -p(n 2 - ul) + pouo{uo -u)-(p~ Po) + P-^ 

2 


(p - po) 


and 


E, = p 


dh 


ds 


( s — So) — Tq{s — so) 


ip 


where 


and 


dh 

ds 


1 kJ 

p r s f 


7-1 


dh 

dp 


-j. i 

= p 1 ST' 


(4.10) 


(4.11) 


(4.12) 


(4.13) 


for an ideal gas. 

It is instructive to examine the energy components and energy source terms in 
space time to see how these quantities change during the sound-shock interaction 
process. These quantities presented in Eqns. 4.10, 4.11 and 4.9 are shown in 
Figures 4.11 -4.13 for the case where the Mach number upstream of the shock, Mi - 
3, the acoustic disturbance amplitude, « = .1. and 512 cells are distributed along 
the duct length. For these calculations, the shock is initially at * = .4. The sound 
disturbance enters the duct at t = 0 at the duct inlet and propagates downstream. 
A, a time of t * 0.13 the sound wave hits the shock. After the interaction, 
Figure 4.11 shows that acoustic energy is present downstream. The amplitude 
of the sound energy downstream of the shock is higher than the energy upstream, 
indicating that acoustic energy is generated by the sound-shock interaction process. 
It is comforting to note that the slope of the path of the sound wave in space-time 
increases after interact, on with the shock. The inverse of the slopes before and 
after the shock corresponds to the quantities », + c, and u 2 + o, respectively, 
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Figure 4.11: Disturbance acoustic energy as a function of space time. Pre-shock 
Mach number is 3, disturbance acoustic amplitude = 0.1, 512 cells distributed 
along duct length. 
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Figure 4.12: Disturbance entropy energy as a function of space time. Pre-shock 

Mach number is 3, disturbance aconstic amplitude = 0.1, 512 cells distributed 
along duct length. 

where the subscripts 1 and 2 refer to the states npstream and downstream of the 
shock. 

Figure 4.12 show, the component of entropy energy in space-time. Clearly, 
because the sound wave is by definition isentropic, them is no disturbance entropy 
upstream of the shock. After the sonnd wave hits the shock at t * 0.13, however 
entropy energy appear, downstream. Thus, entropy energy is generated during 
sound-shock interaction. The inverse slope of the path of the distnrbance entropy 
correspond, to the downstream convection velocity, as expected. Figure 4.13 is 
interesting became it portrays the variation in the source term of E q n. 4.9, and 
provides insight into the nature of the generation of distnrbance aconstic an! en- 
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Figure 4.13: Disturbance energy source as a function of space time. Pre-shock 

Mach number is 3, disturbance acoustic amph.ude = 0.1, 512 cells distributed 
along duct length. 

trcpy energy downstream o, the shock. The source term is aero except along the 
rHock wave for t > 0.13. Thus, the shock wave is the source of the disturbance en- 
erg.es downstream. Inspection of E q ». 4.9 indicate, that the source of disturbance 
energy the transfer of energy from the mean flow. Thi. transfer can occur in the 
presence of a temperature gradient with fluctuations in entropy and momentum 
aud with fluctuation of entropy and temperature. Note that the source term ,s 

aero when the fluctuation in the entropy is aero. Thi. implies that the disturbance 
energy source is related to the shock motion. 
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4.5 Concluding Remarks 

This chapter describes a numerical investigation of sound amplification by a shock 
wave. The MacCormack and 2nd, 3rd, and 4th order ENO schemes are employed 
to compute the time dependent shocked flow field in a converging-diverging noznle. 
The How is disturbed by introducing a sound wave at the inflow boundary. AU of 
the schemes are shown to do a good job in predicting the perturbation amplitude 
and phase speed in the noszle, especially in the supersonic, smooth portion of the 
flow The high-order ENO schemes provide the best overall results because the flow 
around the shock does not contain spurious oscillations and the How downstream 
of the shock shows little dissipation. 

The numerical results are compared to the Unear theory, linear theory vah- 
dates the numerical solution for small perturbation ampht.de. Numerical results 
for the larger perturbation amphtude also compare weU to the Unear theory, indi- 
cating that the Unear theory has a wide range of apphcabiUty in predicting sound 

amplification by a shock. 

Analysis of the equation governing the perturbation energy shows that distur 
bance energy is generated at the shock. The source term for this energy goes to 
ze ro as the entropy fluctuation goes to aero. This indicates the importance of shock 

motion on the generation of disturbance energy. 
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Chapter 5 


Interaction of a Vortex Ring with 
a Shock Wave 

5.1 Introduction 

In this chapter the interaction of a vortex ring with a normal shock wave is consid- 
ered. The study of this interaction is meaningful because it models a fundamental 
mechanism of sound generation in supersonic jets: the interaction of turbulence 
with shock waves. 

Early research on shock vortex interaction focused primarily on experimental 
studies ([1], [2], [3]) and the development of predictive linear theories ([4], [5], [6] ) 
that were compared with experimental results. Ribner suggested that the study of 
the interaction of a vorticity wave with a shock would provide useful information 
regarding the sound generated by turbulence in supersonic jet flows [4]. Ribner 
studied the problem analytically and developed a linear theory which described 
sound generated by vortex-shock interaction ([4], [6]). At about the same time, 
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Moo. studied the factions of a variety of plane wave disturbances with shocks 
in an unsteady reference frame [5], 

Pao and Salas were the first to study the interaction of a shock wave and a vortex 
numerically [7). This study investigated the interaction of a columnar vortex with 
a normal shock wave by solving the Euler equations using MacCormack's scheme 

W ‘ lh “ ShOCW “ mg " UmCTCal leCh " i "“- “»■ Zang, and Hussaini [8], Hussain., 
et al. [9j, and Kopriva et. al [10] applied spectral methods with shock fitting 

accuracy of the solution, but were limited to weak shock-vortex interaction cases 
Meadows, Kumar and Hussaini (11) studied the interaction of a columnar vortex 
with a shock wave using a shock capturing scheme. Shock-capturing proved to 
be beneficial because strong shock-vortex interaction cases which result in the 
formation of secondary shocks, could be studied readily. The authors noted that 
m order to provide a quantitatively accurate representation of the acoustic wave 
improved downstream boundary conditions and higher-order numerical schemes 
were required. Casper [12] then investigated the shock-vortex interaction problem 
with a high-order ENO scheme and found that higher accuracy greatly improved 
the resolution of the acoustic wave downstream of the shock. 

The work cited above has been for columnar vortex-shock interactions. This 
^dy is believed to be the firs, for the interaction o, a vortex ring with a shock 
wave. The interaction of a vortex ring with a shock more closely models the 
interaction of turbulence within the shear layer with shock waves present in the 
Plumes of imperfectly expanded axisymmetric snpersonic jets. 
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Figure 5.1: Vortex Ring - Shock Interaction. 

5.2 Model 

5.2.1 Overview 

The calculations presented in this section are designed to model the interaction of 
a vortical structure within the supersonic portion of the shear layer with a shock 
wave. The core size is small relative to the size of the ring in these calculations, 
primarily because of numerical considerations. Because the shear layer in a jet 
spreads with increasing distance from the nozzle exit plane, a vortex ring with 
a small core models the interaction of disturbances in the shock cell closest to 
the nozzle lip. When results are presented in dimensional units, the variables are 
dimensionalized with ambient stagnation sound speed and atmospheric pressure; 
thus, the calculations performed most closely model disturbances interacting with 
the first shock cell of an overexpanded nozzle, illustrated in Figure 5.1. 

The strength of the vortex is chosen to correspond to the observed strengths of 
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turbulent fluctuations in the jet shear layer. Experimental observations indicate 
that an appropriate level of velocity disturbance is on the order of three percent 
of the mean flow velocity [22]. The vortex strength chosen for the majority of the 
calculations presented here (Mi = 1.5) has a velocity perturbation (defined as the 
ratio of the vortex core velocity to the upstream mean flow velocity) of 6.8%. 

Although imperfectly expanded jet flows typically contain systems of oblique 
shocks, this model uses the interaction of a ring vortex with a normal shock wave. 
This is a reasonable approximation because the oblique shocks tend to curve up- 
stream with increasing distance from the jet centerline and are approximately nor- 
mal upon termination in the shear layer [22]. 

The sense of the vortex rotation is taken to be counterclockwise for the majority 
of the calculations because the velocity decreases with increasing distance from the 
jet axis, resulting in vorticity of a positive sense. Some results are shown for the 
case of a clockwise rotating vortex in Section 5.5. 

5.2.2 Geometry 

The geometric model of the interaction of a ring vortex with a shock wave used in 
the computations is illustrated in Figure 5.2. A vortex ring is introduced upstream 
of the shock wave. The vortex ring is characterized by its strength, T, its core 
radius, r c , and the distance from the axis of symmetry to the vortex filament, tq. 
The shock is characterized by the Mach number of the flow upstream of the shock, 

Mi. 
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F {/>n,pu ! + p t p UVi + 


(5.3) 


G-{pv, pm, ,,v‘ + p _ ( ?e + p)vj T 


(5.4) 
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H = {0,0,p,0} t (5.5) 

and standard notation is used for the flow variables. Thus, p is density, u and v 
are the velocity components in the axial and radial directions, e is total energy per 
unit mass, p is pressure, x is the axial position, and r is the radial distance from 
the axis of symmetry. 

5.2.4 Boundary Conditions 

Conditions are prescribed at the inflow and outflow computational boundaries to 
establish a shock with the specified strength at x = 7. The boundary condition 
along the axis of symmetry is prescribed by requiring that the fluxes across the 
axis are zero. The method of prescribing the boundary conditions is described in 
[13]. 

5.2.5 Solution Procedure 

The calculation is performed in two steps. First, a steady, shocked flow is es- 
tablished with the flow parallel to the axis of symmetry. Next, a vortex ring is 
introduced to the flow field 7 core radii upstream of the shock at T = 0. All flow 
variables are normalized with respect to the static pressure and density upstream 
of the shock, and the size of the vortex core radius. The vortex ring then convects 
downstream, passing through the shock wave. After this interaction, significant 
changes occur downstream of the shock wave. Sound, vorticity, and entropy are 
generated as the vortex interacts with the shock wave. The primary focus of this 
work is on the acoustic waves generated by the interaction. 
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5.2.6 Vortex Model 


A model of the ring vortex is reared to introduce appropriate perturbarions 
rn velocity and pressure to the flowfield as an initial condition. The vortex is 
introduced at the initial time t = 0, and at subsequent times the Euler equations 
determine the entire how held - including the vortex. A, all subsequent time steps 

EUkl e<,Uali0nS *" “ corresponding to the truncation 

error of the ENO algorithm, which is third order accurate in space and time. 

The initial condition prescribed for the vortex is the classical toroidal vortex of 

Lnmb [14], augmented with a sobd body core. The derivation of the velocity and 
pressure equations is provided in Appendix B. 

Limits of the Vortex Ring Model 

The cross section of the vortex ring core is approximately circular when the vortex 
filament is 170 core radii away from the axis of symmetry [16], When the vortex 

core radius increases relative to the ring radius, the degme of circularity decreases, 

until the vortex becomes Hill’s vortex flfil an A a 

vortex [lojj and the streamlines become oblong and 

flatten close to the axis of symmetry. To closely approximate the vortex ring with 

symmetry un.ess otherwise noted. The small siae of the vortex core radius relative 
40 the ring radius reasonably approximates thesiaeof a turbulent structure within 
the supersonic portion of the shear layer near the first shock cell. 

5.2.7 Vortex Parameter Modeling 

The parameters of the vortex ring which model the physical jet turbulence/shock 
interaction will be described here. The first parameter of interest is the vortex 
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strength. The circulation is typically considered to be a measure of the vortex 
strength. For the vortex ring, the circulation is related to the geometrical param- 
eters and induced velocity, V, according to 116] 

_ 4Vjrro (5.6) 

The next parameter of interest is the ratio of the core radius to the ring radius, 
I*. The core radius models the size of the disturbance in the supersonic portion of 
the mixing layer interacting with the shock wave, while the ring radius models the 
radial distance between the jet axis of symmetry and the inner boundary of the 
turbulent mixing layer. Obviously, this ratio is a function of distance downstream 

from the jet exit plane. 

5.2.8 Vortex Preservation Study 

An important feature of this calculation is the accurate representation of the source 
of the sound generation: the interaction of a vortex with the shock. To ensure that 
the vortex is adequately resolved in the calculation, the time history of the mini- 
mum pressure within the vortex is tracked. A series of computations is performed 
where the vortex ring convects over grids of various resolutions. It was found that 
while convecting on a uniform mesh with 10 cells per core diameter, the minimum 
pressure varies by a maximum of 0.035% over the time it takes the vortex to travel 
7 core radii, which is the initial distance of the vortex from the shock. This level 
of numerical error is considered to be acceptable, and so this grid resolution was 


held fixed for all the computations. 
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5.2.9 Computational Grid 

As mentioned in the previous section, it was found that a grid having 10 cells 
per vortex core diameter was sufficient to adequately resolve the vortex in the 
computation. This grid resolution is used in the fine uniform grid in the region 
—6 < x < 85 and 72 < r < 178, where Ax = A r = 0.2. However, it is computa- 
tionally prohibitive to use a uniform mesh with such fine resolution throughout the 
entire computational domain. The computational domain is required to be large 
enough such that at least one wavelength of the acoustic disturbance is contained 
in the calculation. This allows accurate measurements to be made in the acoustic 
far field so that computations to determine the acoustic energy level can be made. 
In addition, the computational domain must also be large enough that once the 
acoustic wave has passed a far field observer point, the wave has time to go through 
one period of its oscillation without being contaminated by errors introduced by 
the boundaries. Thus, in the calculations shown here, the computational grid is 
uniform over the range of the acoustic source interaction and sound wave prop- 
agation through at least one acoustic wavelength. The grid is then stretched to 
the boundaries using a hyberbolic sine function. The grid contains 481 cells in the 
axial direction and 556 cells in the radial direction. Every 15th cell of the grid is 
shown in Figure 5.3. 

It is important to note that this grid resolves the generated acoustic waves. 
Preliminary test calculations performed over Mach number range 1.1 < M < 
1.7 cover an acoustic peak-to-valley wavelength range of 2.5 < A < 3.0. Thus, 
the minimum resolution of the acoustic wave is 12.5 points per peak-to-valley 
wavelength. The actual number of points per wavelength is at least twice this 
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Figure 5.3: Standard grid used in calculations. Only every 15th cell is shown. 

number. Thus, the acoustic wave is well resolved. 

It is also important to note that there is another length scale of possible interest 
in this problem, which is not well resolved on this grid. During the interaction of 
the vortex with the shock, the shock wave begins to move. For typically small vor- 
tex strengths, the shock displacement is small. In fact, for many of the calculations 
performed during the course of this research, the shock moves through the distance 
of only one cell. Because further grid refinement proved to be prohibitively expen- 
sive and one-dimensional calculations of sound-shock interaction provided results 
which matched linear theory well without resolving the shock motion, the grid was 
not refined further. 
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5.3 Typical Interaction- Counterclockwise 
Vortex 

5.3.1 Overview 

A typical interaction of a counter-clockwise (CCW) rotating vortex with a shock 
wave is is described in this section. In the calculations for this case, the Mach 
number upstream of the shock is M = 1.5, and the vortex strength is T = .75. 
In the contour plots shown in this chapter, contours are shown for a part of the 
computational domain near the interaction point of the vortex core and the shock 
wave. Because of the axisymmetry, only a single cross-section of the solution will 
be shown in the results which follow. All the results presented are within the 
uniform grid region of the computational grid so that effects of grid stretching are 
not present in the results. The contours are computed on a grid which is five times 
as coarse as the computational grid. Note that the contour range is kept constant 
for each plot so that relative values of the perturbation quantities at different time 
levels can be compared. In these calculations, a unit of time is defined as T = r c ju\ 
where u\ is the upstream flow speed; thus T represents the time it takes the vortex 
core to move 1 radius in the flow upstream of the shock. 

Contour plots of the flow variables are shown for three selected times: T = 0, 
T — 8, and T = 50. At time T = 0, the vortex is upstream of the shock. At 
T = 8 the most upstream edge of the vortex core is approximately aligned with 
the shock; and at T = 50, the vortex is approximately 30 core radii downstream 
of the shock. 
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5.3.2 Pressure 

Figure 5.4 shows the change in pressure from the mean state. The only perturb 
hr pressure a, T = 0 ,s the decrease pressure a, the vortex core. As the vortex 
begins to interact with the shock, additronal pressure perturbatrous are generate 
downstream of the shock. At T - S Mgh amphtude pressure disturbances are seen 
JU st downstream of the shock. As time passes, these pressure drsturbauces travei 
downstream more rapidly than the vortex. Figure 5.5 shows a plot of the pos.tron 
of the peak pressure (in radri from the vortex filament) a, a functron of t.me (m 
periods). The plot shows a linearly increasing change in the position of the pressure 
disturbance as a function of time, which means that the pressure disturbance rs 
traveling at a constant speed relative to the vortex. The slope of tins curve ,s foun 
to be Af = 0.77 core radii per period, the sound speed downstream of the shoe . 
(The velocity 0.77 core radri per period is equivalent to a velocity normalrsed by 
the sound speed upstream of the shock of U5. The ratio of the sound spee s 
across the shock is 1.15 for an upstream Mach number of 1.5. ) Thus, the pressure 
disturbances are traveling a, the sound velocity relative to the mean how, satrsfymg 

one of the defining features of acoustic waves. 

The structure of the acoustic wave is readily apparent at T = 50 m Figure 
5 4 As predicted in linear theory, the sound wave is guadrupole in nature: the 

• d alternate compression-rarefaction-compression- 

acoustic wavefront is comprised of alternate comp 
rarefaction fronts. 

Pressure perturbations along radii extending from the vortex center at 10 degree 

increments are shown in Figure 5.6. This figure quantifies the change m amphtude 
as a function of the angle from the horizontal, and shows that the peak-to-pea 


108 


ORIGINAL PAGE 

Cl L Oft PhC TC'C r: APf-j 



Figure 5.4: Contours of 
T = 8, and T = 50. 


pressure perturbation downstream of the shock at T = 


0 , 


pressure amplitude i, marimum at # = 50 and # = -55 degrees. fVom this figure, 

the valley.to.peak measure of the wavelength of the pressure disturbance, is , Mn 
to be A ta 2.8. 

Figure 5.7 show, a carpet plot of the pressure levels downstream of the shock 

* r Th ' “‘ 8 ' " Pr0 “’ Md Wilh »■* [24] and is shown at an angle to 

danfy the detailed feature, of the fiow field. Note the ^solution of the cylindrical 

the cylindrical acoustic wave front. 

Figure 5.8 show, the product of the , q uare root of the distance the acoustic 
wave ha, traveled (relative to the verier core) and the peak pressure amplitude 
along a line extending from the verier core a, 60 degree, from the horizontal 
passing through the verier fiUmen, a. a function of distance from the verier 
filament. The product varies significantly a, small distance, from the core, which 
.. anonymous with early time and prorimity to the source, but flatten, to an 
aw constant level at larger distances. This show, tha, in the far field of these 


109 


Distance 

from 

core 



Figure 5.5: Position of the peak pressure perturbation (in core radii from the vortex 
filament) as a function of time (in periods). The slope of the curve is the sound 
speed downstream of the shock. 

calculations, the sound pressure decays as which is characteristic of geometrical 
spreading in two dimensions, commonly known as cylindrical spreading. Although 
these calculations are performed in an axisymmetric coordinate system which, in 
general, should allow for three dimensional geometrical spreading of the sound, the 
presence of the shock wave prohibits the spreading upstream. In realistic flows, 
where viscosity plays a role in the fluid dynamics, shock waves are finite in extent. 
In these flows, cylindrical spreading would occur close to the shock wave, while 
spherical spreading (with p decaying as £) would occur in the far field. Hardin [23] 
has shown that in cylindrical coordinates, the sound decay rate is a function of the 
source size to the distance from the source to the observer. Later it will be shown 
that the interaction of the vortex ring and shock wave produces a disturbance 
which travels along the shock, thereby increasing the size of the potential noise 
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Figure 5.7: Pressure perturbations downstream of the shock at T - 50 
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constant value, indicating that the 


pressure magnitude is shown to 


asymptote to a 


acoustic wave spreads cylindrical^. 
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seven core radii upstream of the shock for most of the calculations presented in this 
paper. However, to test the effect of the initial distance between the vortex and the 
shock wave, A*i, on the flow field downstream of the shock, a single calculation 
is presented where the vortex is initially twenty-three core radii upstream of the 
shock. Figure 5.9 shows the solution downstream of the shock at times T = 50 
and T = 66 for initial shock locations of x = 7 and x = 23, respectively. In both of 
these solutions the vortex has traveled approximately 23 core radii after its filament 
has passed through the shock wave. The cylindrical portion of the pressure field is 
similar for both initial conditions; however, there are differences in the downstream 
flow-field. Most notably, the cusp-shaped pressure structure is much larger for the 
longer time solution (A*i = 23), indicating that it originates at the beginning of 
the computation. Further study of this feature is required to quantify the effect of 
axisymmetry and the initial condition on its structure, and to determine whether 
it is numerical or physical in nature. 

5.3.3 Density 

Figure 5.11 shows the change in density of the flow field. At time T = 0 the 
only change in density is associated with the vortex. As the vortex core begins 
to interact with the shock, density waves appear downstream of the shock. The 
interesting character of the density waves becomes clear in the final figure of the 
time sequence. Two types of density perturbations are evident in the figure. Den- 
sity disturbances associated with the acoustic wave propagate in a nearly circular 
pattern downstream of the shock. This is to be expected, because acoustic waves 
are by definition isentropic and there is a clear relationship between the density 
and pressure (s = constant — In addition, there are also density disturbances 
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' gUre 5 ' 9: PreStUrc COntOUrE d ° Wnstrea “ for a time T = 50 (left figure) and 

11 = 66 (rigi ‘ ^ iE >«*- - approximately twenty., hree core 

radn down,, ream of the shock for both cate. The difference in the relation, i. a 

re.nl, of the initial placement of the vortex relative to the .hoch. In the figure on 

the left, the vortex i. initiafiy .even core radii away from the .hoch. In the ff g ure 

on the right, the vortex is initially twenty three core radii away from the shock. 
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of a. convective nature. These density disturbances look like spokes reaching out 
from the vortex core and terminating at the shock. These disturbances are asso- 
ciated with the entropy waves associated with changes in the shock strength. As 
the vortex interacts with the shock, the shock wave begins to move. As the shock 
moves, the change in entropy across the shock is no longer constant, and an en- 
tropy wave is generated which convects downstream at the local flow velocity. For 
the counter clockwise rotating vortex, the portion of the shock above the vortex 
filament initially moves upstream and the portion of the shock below the filament 
moves downstream. As time increases, disturbances move away from the interac- 
tion point, and their motion along the shock wave creates entropy perturbations 
which convect downstream at the local flow velocity. These disturbances show a 
strong resemblance to the features visualized by Naumann and Hermanns [3] in 
their experiment and described as contact surfaces. 

Figure 5.10 shows a digitized version of a Mach-Zehnder interferogram from [3] 
which illustrates the sound wave and contact surface visualized in the experiment. 
In the experiment, the interaction is produced as follows. A sharp edged profile is 
placed in a shock tube. A diaphragm (located to the left of the airfoil) is broken and 
a weak disturbance travels through the tube and over the trading edge, producing 
a starting vortex. Another diaphragm (located at the right of the airfoil ) is 
broken and a shock wave travels towards the vortex. The results of the interaction 
are shown in Figure 5.10. Because the interaction observed in this experiment is 
strong, the motion of the shock wave is pronounced and the acoustic wavefront is 
asymmetrical. The contact surfaces form a curved funnel-shaped structure between 


the shock wave and the vortex. 
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Figure 5.10: Mach-Zehnder interferogram of sound wave and contact surfaces gen- 
erated shock- vortex interaction. (From [3]) 
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Figure 5.11: Contour, of denrity perturbation down.tream of the .hock et 
T = 0,r = 8, and T = 50. 

Figure 5.12 .how. . carpet plot of the den.ity at T = 50. The .ound wave and 
complex nature of the contact .urfhce. me clemly virible. 

Experiment, provide evidence for the phyrical nature of the contact .urface. 
observed in the.e computation.. However, earher analyri. in Chapter 3 demon- 
■trated that the computation of .lowly moving .hock wave, produce, error which 
manifest* it. elf primarily in entropy. Therefore, caution mu,t be exercried in the 
interpretation of the strength of the.e di.turbance.. Further analyri. i. required 
before the nature of these contact surfaces is validated. 


5.8.4 Vorticity 


Figure 5.13 .how. the change in vorticity. A. the initial time, the only vorticity i. 
• cwcular .pike at the vortex core, so this contour plot i. no, shown. A. the vortex 
interacts with th, .hock a, T = 8, vorticity appear, down., ream. A, T = 50, 
the character of the vorticity become, clearer. The vorticity pattern., like the 
convective density dirturbance, look like the .poke, of a wheel radiating horn 
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Figure 5.12: Contours of density perturbation downstream of the shock at T 


50. 
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Figure 5.13: Contours of vorticity perturbation downstream of the shock at T = 8 
and T = 50. 

the vortex core and terminating at the shock. Zooming in on a region in the 
immediate vicinity of the vortex core at T = 50 (Figure 5.14), it is clear that the 
lines of constant vorticity are oblong in shape. 

5.3.5 Velocity 

Figures 5.15 and 5.17 show the changes in the axial and radial components of veloc- 
ity. Except for small changes in amplitude, the velocity field shows no significant 
variation as a result of the interaction on this scale. There are, however, very in- 
teresting features similar to those observed in the density perturbations when the 
contour range is decreased. Figure 5.16 shows the change in the axial component 
of velocity on a smaller contour range. The vortex, acoustic wave, and entropy 
disturbances are all readily distinguished on this scale. 

5.3.6 Entropy 

Figure 5.18 shows the change in entropy, defined as 8s — s a ^ — sq. At the 
initial time, there is no fluctuation in entropy, since the initial vortex is isentropic. 
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Figure 5.14: Contours of vorticity perturbation 
rounding tie vortex filament at T = 50. 
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Figure 5.16: Contours of axial velocity perturbation at T — 50. The range of the 
contour levels has been reduced to show the interesting velocity features. 



Figure 5.17: Contours of radial velocity perturbation downstream of the shock at 
T = 0, T = 8, and T = 50. 



122 


r 


160 

140 

120 

100 

80 



x 




Figure 5.18.- Contours of entropy perturbation at T = 8 and and T = 50. 
Therefore, the contour plot at 7 1 = 0 is not shown. However, as tie vortex core 

“ eraC ‘ S mti ,hC Si ° Ck “ tr ° W a PP eara downs, reanr. Note that the entropy 

perturbations correspond to the convective fluctuations in density and vorticity as 
seen in Figures 5.11 and 5.13. 

5.3.7 Discussion 

Figures 5, through 5.18 illustrate interesting How features which result from the 
mteractum of a nng vortex and a shock wave. After observing the perturbation. 

" “ ‘ 8 *“ ** “ ~ -e and contact surfaces result 

, tUS m ‘“' The “ "" ** P-P-y that it is isentropic, and 

-els at the sound speed relative to the moving fluid. The contact surfaces do no, 

support a pressure disturbance but do support differences in flow velocity, density 

entropy, and vorticity. The disturbances are the result „f ,1. • . 

tte result °f ‘he interaction. What is 

' ‘ OU ‘ in ‘ eraCti0n WUdl Crci,eS disturbances? Prom a mathematical 
perspective, the shock wave is the feature which allows for the coupling of the 

near mode, (acoustic, entropy, and vorticity). Crocco, Theorem provides the 

relationship between the flow vorticity and entropy, and from this theorem it can be 




123 


200 

150 

100 

50 




6X. 

0.16 

0.13 

0.10 

0.07 

0.04 

0.01 

• 0.02 

- 0.05 

- 0.07 

- 0.10 


Figure 5.19: Shock displacement as a function of space-t im e 

shown that vorticity can be generated in the presence of a curved shock. Because 
the evidence points to the importance of the shock wave in the generation of 
the acoustic, vortical, and entropy disturbances, a closer observation of the shock 
dynamics will be made in the next section. 

5.3.8 Shock Dynamics 

Figure 5.19 shows a space time diagram of the shock displacement. For times 
T — 0 through T = 6 the shock position remains nearly fixed as the vortex convects 
supersonically towards it. Once the vortex core begins to strongly interact with 
the shock at T = 6, the shock motion becomes significant. Figure 5.19 illustrates 
the nature of the shock displacement as a function of time. A small disturbance on 
the shock is initiated at T ~ 6. The disturbance proceeds to travel both away from 
and towards the axis of symmetry. The speeds of these waves is approximately the 
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e «. of symmetry, be, continues to grow towards the ax, of symmetry. 

F.g„re 5.20 is included to clarify the detaiis of the shock displacement at various 

t,meS ' fi8Ur ' SWS the ““ d “ Pl — -* - ** shock as a function of radial 
P°s,t.o„ along the shock for T = 1, r , 6 , r = 10 , r _ 20> T ^ ^ T = ^ 

^0 I ~ 50 MT= 10 ’ ^ 8hOCk ^ ^ ma3dmUm ^Placement magnitudes 
° “ 014 “ the Upslream and downstream directions, respectively. At this 

ear y tune, the displacement of the shock is limited a small region around r - 

125 the radial position o, the vortex filament. As time progresses, the maximum 

SPkCeme "‘ ° CC “ rS interaction point. The maximum 

upstream displacement occnrs a, a time of T Ri 10. However, the displacement 

“ 6 t '°” Stream direCli °” ™ • even after 50 periods, ff one 

considers the motion of the maximum shock displacements along the shock wave it 

,S ” lhat lheM from the vortex-shock interacts 

Pomt at a particular velocity. By careful examination of the results presented 

" Figure 5.10, this velocity is found to be « 0.5d core radii per period (which 
corresponds to a normalised velocity o, « This speed is well below both the 
upstream and downstream sonnd speeds (1.0 and 1.15, respectively, 

I« « interesting that the shock wave displacement continues to increase in the 

direction towards the center of . 

of the ring, while it maximises and then decreases in 

the direction away from the ring’s axis of symmetry. 

Figures 5.10 and 5,0 show the asymmetry in the shock position relative to 
an axis passing through the vortex filament position a, r = 125. The portion of 
t e shock closest to the axis of symmetry continues to move downstream over the 
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tune of the calculation. The maximum shock displacement is greater downstream 
(A®, _ 0.16) than upstream (Ax, = -0.1). 

Figure 5.21 shows the displacement of the shock and the perturbation in density 
and pressure immediately downstream of the shock (a. a = 7.02) at T = 50. This 
figure shows that there is a correlation between the disturbances downstream of the 
shock and the shock displacement. The shock displacement leads the disturbances 

“ “ “ d ~ <>f ‘he shock displacement 

correspond to the large Jumps in density at r = 100 and , = 150. The changes in 

slope of the fix. vs , curve located at , = 70, r = 120, and r = 130 also correspond 

to significant high frequency changes in density. This type of fluctuation in pressure 

occurs only a, r = 70. Thus, by comparing this figure with the contour plots of 

Figures 5.4 and 5.11, the spoke-like disturbances in density and vorticity can be 
related directly to the shock motion and curvature. 

5.3.9 Frequency Analysis 

In order to obtain information regarding the frequency content of the solution, the 
Fourier transform is used to transform the time dependent information obtained 
from the numerical calculation to the frequency domain: 


/ -f-oo 

Sp(t)e- iut 

-oo 


dt 


where fip is the perturbation pressure defined as the difference between the 
and the mean pressure (fip = p - *), and „ is the cyclic frequency. 

This integral is written in discrete form: 


(5,7) 


pressure 
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Figure 5.21: Shock displacement and density and pressure perturbations as func- 
tions of radial distance (on the vertical axis) for T = 50. The density and pressure 
perturbations are obtained slightly downstream of the shock (* = 7.02). 
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N - 1 

PxAtJ2 p(nAt)e~ iunAt (5.8) 

n=0 

where At is the discrete time interval of the discrete pressure time history, p(nAt), 
N is the number of points in the time history. The summation Equation 5.8 is 
computed using the Fast Fourier Transform (FFT) approach. In order to efficiently 
use the FFT, N must be a power of 2. 

In this numerical simulation, solutions were saved at every period (defined 
previously as T = ^-). Thus, there are 50 discrete representations of pressure in 
the time history available for analysis at each point in the computational domain. 
In order to effectively use the FFT approach, the time history at each point at 
which the spectrum is to be computed is padded with zeros so that the number 
of discrete values in the time history is a power of two. Padding the time history 
with zeros does not affect the highest frequency which can be resolved in a discrete 
approximation to the Fourier transform (the so-called Nyquist frequency), since 
the Nyquist frequency u> c is a function only of the temporal increment in the time 
history, u> c = However, padding the time history with zeros does increase the 
frequency resolution of the spectral estimate. 

The increase in frequency resolution provided by zero padding is readily un- 
derstood by an example [21], Suppose the original time history T has b elements: 
T(nAf), n = 0,1,2, ..b — 1. The original time history is then padded with an 
additional b elements, which are zero: T(nAt) = 0, n = b,b + 1, ...26 — 1. The 
discrete Fourier transform Ty is then 

TtM = £ 2b E T{nAt)e~ i2Tkn l 2h 

71=0 

A . 26—1 . , ,, 

= ^ E T(nAt)e~ %Tkn l b 

n—0 


(5.9) 
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Note that w* = k = 0,1,2, ...b, and the frequency resolution is Au; = 5^. 

Similarly, for the original time history without zero padding, w* = k — 

0, 1, 2, ...6/2, and the frequency resolution is Aw = Thus, the zero padding 
has the effect of increasing the frequency resolution by a factor of two. 

Although it would be sufficient to pad the 50 element time history with 14 zeros 
to obtain a 64 element time history, each time series was padded with 206 zeros for 
a total of 256 points in the time history. This provided better frequency resolution 
necessary for localization of the acoustic energy in the spectrum. The price paid 
for enhanced frequency resolution is a loss of accuracy. The reader is referred to 
[21] for details. 

Typical sound pressure levels are presented in Figures 5.22 and 5.23. The 
frequencies on the horizontal axis are normalized by the upstream flow velocity 
and core radius, to obtain a Strouhal number. In these figures, the effect of the 
mean flow has been removed from the spectra by subtracting the time average of 
the pressure from each point in the time history. This removes the energy from 
the zero frequency bin. The sound pressure level is represented as decibels (dB), 
per convention. The decibel level is obtained by computing: 

P(u)[dB] = 20 log(^^-) (5.10) 

Pref 

where p re f = 20 fiPa is the conventional acoustic reference pressure. 

Figures 5.22 through 5.24 show the sound pressure level at locations along a 
45 degree line from the sound source point (x = 7, y = 125) at 6.08 3 11.74 and 
23.1 core radii away from the source point. The highest energy levels are in the 
low frequency range (Strouhal numbers less than 0.1). There are additional peaks 
at Strouhal numbers of « 0.2, 0.4, 0.6. The energy in the frequency band of about 
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0.22 becomes relatively more important than the low frequency energy as distance 
from the source point is increased. By observation of the time history from which 
the spectrum is computed, the period of the acoustic wave is determined to be 
approximately 21 periods. This corresponds to a Strouhal number Of U;Tr - ZTTfr 

2r(l ) Ul ~ Tui — 

(21X1-5) = °- 2 ' Thus ’ the ener gy in the frequency range of 0.22 is associated with 
the acoustic energy of the signal. The large peaks near the Nyquist frequency 
(Strouhal No. = tt ) are probably due to aliasing. Hardin [20] has shown that 
for spectral estimates near the Nyquist frequency, u> c , contributions from -u c can 
appear, even for a sufficiently sampled time history. Figure 5.25 shows the sound 
pressure level of the two frequencies 0.025 and 0.2 as a function of distance from the 
sound source point. The data are taken along a 45 degree line from the horizontal 
passing through the vortex core at r = 125. The data shows a rapid decay in the 
energy associated with the low frequency energy. The energy decays at about MB 
per doubling of distance away from the source. This implies a ^ decay rate, which 
is consistent with the expected decay rate for a vortical pressure field. 

The energy at the Strouhal number of 0.22 decays at a lower rate % MB per 
doubling of distance, which implies a sound pressure level decay of A. Because 
the sound pressure level is a function of the square of the pressure, this result 
is consistent with the pressure decay rate of ^ observed in Section 5.3.2., and 
implies cylindrical spreading of the acoustic energy. 

5.3.10 Sound Intensity Level 

The sound intensity is a vector quantity defined by: 

T(D= l -£wit 


(5.11) 
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Figure 5.22: Sound pressure level (SPL) pressure as a function of dimensionless 
frequency. Distance from the source is 6.0811 core radii. Angle from horizontal is 
45 degrees. 
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Figure 5.23: Sound pressure level (SPL) as a function of dimensionless frequency. 
Distance from the source is 11.74 core radii. Angle from horizontal is 45 degrees. 



132 


SPL 

(dB) 



Figure 5.24: Sound pressure level (SPL) as a function of dimensionless frequency. 
Distance from the source is 23.1 core radii. Angle from horizontal is 45 degrees. 



Figure 5.25: Sound pressure level (SPL) as a function of distance from the inter- 
action point. This data is taken along a line at 45 degrees from the point where 
a horizontal line passing through the vortex filament (r=125) passes through the 


undisturbed shock. 
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where T is the period of a periodic disturbance or a sufficiently long time for 
nonperiodic bounded signals, and W is the instantaneous acoustic energy flux 
vector: 

W — ( Sp + poUo • Su) 

where 6p, 8u, and 6p are disturbances in pressure, velocity, and density, defined: 

8p = p — p o 
6u = u — u o 
8p = p — p o 

and the subscript 0 represents the mean flows state. 

Note that in a quiescent medium (where Uq = 0), three of the terms in the in- 
stantaneous acoustic energy flux vector vanish, and the equation for sound intensity 
reduces to: 

/.(*) = ^ J*8p8udt (5.12) 

Because the early theoretical development of acoustics was developed for sound 
disturbances in a quiescent state, this definition will be called the “classical” defi- 
nition of sound intensity, and it will be labeled with a subscript c to distinguish it 
from the definition of intensity presented in Equation 5.11. It will be shown later 
that it can be beneficial to consider the classical definition of sound intensity, even 
in the presence of mean flow. 

In the calculations performed here, the acoustic signed is transient, which pro- 
vides some difficulty in the interpretation of sound intensity. Because this is a 
problem which models the periodic convection of turbulent disturbances through 
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is not sound intensity, since most of the disturbances related to the vortex convect. 
The high intensity level shown in this region demonstrates the high correlation 
between the pressure and velocity disturbances in the vortex. (Some acoustical 
energy may be generated as the vortex changes shape after its interaction with the 
shock, but small scale vortical motions were not observed to be a significant sound 
source in these calculations. This is most clearly seen in Figure 5.7 which shows a 
carpet plot of pressure fluctuations downstream of the shock. The most significant 
structure is the ring which has already been identified with the interaction of vortex 
core with the shock wave. Much less significant pressure disturbances are seen in 
the center of this acoustic wave, but there are no waves visibly emanating from the 
vortex.) 

The other region of high intensity is seen along the shock wave. The intensity 
plots show that sound is directed primarily in directions closely aligned with the 
shock wave, originating at the point of interaction and traveling both towards and 
away from the axis of symmetry. Both axial and radial components of intensity 
show high levels along the shock; the region is much narrower in the contour plot 
of the axial component. The high amplitude of these waves makes it difficult to 
classify them as acoustic in the classical sense, but because there is no convective 
velocity in the radial direction, these disturbances are clearly not convective. 

A plot of the axial component of sound intensity using the classical definition 
is shown in Figure 5.27. (Because the mean flow is in the axial direction only, 
the radial component of classical sound intensity is the same as that shown in 
Figure 5.26.) The region of high intensity in the vicinity of r = 125 is due to 
the passage of the vortex along this path. Figure 5.27 shows that the intensity 
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Figure 5.26. Sound Intensity Level. I T is radial component. I x is axial component. 


has four lobes not related to the convection of the vortex. Two of these lobes are 
along the shock wave, similar to the results shown in Figure 5.26. However, two 
additional lobes are seen in this figure. These lobes originate at the point of vortex 
filament-shock interaction and point at angles of approximately 50 degrees and -55 
degrees. The lobe directed at -55 degrees is narrower than the lobe directed at 50 
degrees. It is interesting that these regions are not evident in the contours of the 
full description of the sound intensity (Figure 5.26), because the energy associated 
with the acoustics is overwhelmed by the contributions from the mean flow terms. 
The high intensity level along the lobes is significant because it corresponds to the 
location of the strongest region of the sound wave. (See Figure 5.6.) This indicates 
that it can be beneficial to consider the classical definition of sound intensity, even 
in the presence of a mean flow. It is legitimate to apply Equation 5.12 as long as 
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Figure 5.27: Sound Intensity Level I is • i 

cx the axial component of sound intensity 


level 


-tag the classical definition of sound intensity. Note that in this definition 
the primary directivity of the sound wave is along the shock 
at angles of approximately ±50 degrees. 


wave and downstream 
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it is not used to draw conclusions about acoustic energy conservation. 

5.3.11 Effect of Mach Number on Directivity 

In order to determine the effect of Mach number on the directivity of the pressure 
disturbances, a series of computations was performed where all flow parameters 
were held fixed except the upstream Mach number. The pre-shock Mach numbers 
studied were: M = 1.1, 1.3, 1.4, 1.5, and 1.7. These Mach numbers were chosen 
because they are within the range of practical interest. The sound intensity level 
was computed over this Mach number range. The directivity angles as a function 
of shock strength are plotted in Figure 5.28. The figure shows the angles of the 
four primary beaming directions. The change in lobe direction is most significant 
for the lobes closely aligned along the shock. As the Mach number increases, these 
lobe angles approach ±90 degrees. 

5.3.12 Effect of Flow Mach Number on Sound Pressure 
Level 

Because of the complexities in defining and computing the sound intensity level 
for a transient signal in the presence of a mean flow, and because the human ear 
responds to pressure fluctuations, a study of the effect of flow Mach number on 
the sound pressure level is presented here. In this study, sound pressure level is 
computed for a point in the acoustic far field (defined to be one wavelength away 
from the source) at a distance of 23 core radii from the interaction point at an angle 
of 0 = 50 degrees. In these calculations, the sound pressure level is computed by: 


SPL = 20 Log 


Prma 
Pref _ 


(5.13) 
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Figure 5.28: Directivity Angles as a Function of Upstream Mach Number. The 
upper curve corresponds to the intensity lobe closest to the shock above the vortex 
filament; the second curve from the top corresponds to the intensity lobe to the far 
nght of the shock and above the filament; the third line from the top corresponds 
to the intensity lobe to the far right of the shock and below the vortex filament 
position, and the bottom curve corresponds to the intensity lobe closest to the 
shock and below the vortex filament. 
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where p Tm , = sqrt± f 0 T pdt, and where T is the period of the acoustic signal. Be- 
cause this calculation is performed in the far field, the acoustic wave is readily dis- 
tinguished from hydrodynamic disturbances. The result of this study is presented 
in Figure 5.29 which shows sound pressure level as a function of /3 = sjM 2 - 1, for 
Mach numbers in the range 1.1 < M < 1.5. In this figure, the solid fine with small 
dots represents the sound pressure level obtained from the time histones provided 
by the ENO computations of the ring vortex - shock interaction. The large circles 
are the experimental measurements obtained by Seiner and Norum [19] for shock 
noise of an underexpanded supersonic jet, measured at a distance of 12 feet from 
the jet centerline (approximately 146 jet radii from the source). The dashed line 
shows the trend first observed by Harper-Bourne and Fisher [18] that the sound 
pressure level of shock noise in imperfectly expanded supersonic jets is proportional 
to /? 4 over a large range of flow Mach numbers. The slope of the sound presure 
level versus 0 curve very nearly matches the slope of the experimental results for 
shock noise in supersonic jets. This is a significant result because it shows that this 
simple model for shock noise generation can predict the effect of Mach number on 
shock noise observed in experiment. 

Although this model correctly reproduces the trend in sound intensity level as 
a function of shock strength, it does not reproduce the actual sound amplitude, 
even when differences in the distance between the sound source and measurement 
location are accounted for. This is not surprising since shock noise measured during 
an experiment is the result of many interactions of turbulent structures of a variety 
of sizes and strengths with a sequence of shock waves of decaying strength, and at 
varying angles to the flow; this calculation is the result of only a single interaction. 



141 



P 


Figure 5.29: Sound pressure level (SPL) as a function ( 3 , where (3 is a measure of 
upstream Mach number. 
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5.4 Strong Interaction 

The results presented in this section are for the interaction of a strong vortex with 
a shock wave. For all the results presented in this section, the upstream Mach 
number is 1.5 and the strength of the vortex is T = 5.5. Thus, this vortex is l\ 
times stronger than the vortex studied in the previous sections. When the vortex 
strength increases, the pressure disturbances generated downstream of the shock 
increase in magnitude, as illustrated by Figure 5.30 which shows contour plots of 
the pressure disturbances at T = 50. The actual range of the pressure disturbances 
are on the order of 1, but, for the plot, the range of contours was chosen to be 
from -.01 to .007 so that direct comparison with Figure 5.4 can be made. It is 
interesting to compare the results of Figure 5.30 with Figure 5.4. In both figures, 
a sound wave and a cusp-like pressure disturbance are apparent downstream of 
the shock. However, the sound wave is much stronger in the strong vortex-shock 
interaction case (Figure 5.30), and on this plot contour scale, additional pressure 
structures are visible downstream and below the vortex filament . 

Figure 5.31 shows pressure perturbation along radii extending from the vortex 
center at 10 degree increments. This figure shows that the peak-to-peak pressure 
amplitude is a maximum at 55 degrees and -55 degrees, as in the case for T = 0.75. 
From this figure, the valley-to-peak measure of the wavelength is found to be « 2.5. 
By comparison with Figure 5.6, it is clear that the disturbances for the strong 
vortex interaction case are almost an order of magnitude larger. For example, the 
amplitude of the peak disturbance at 50 degrees is 10.2 times larger for the strong 
interaction case (T = 5.5) than the weak interaction case (r = .75). Thus, the 
acoustic pressure p scales as the vortex strength T. 
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Figure 5.30: Contours of pressure perturbation downstream of the shock at T = 50. 
Vortex rotation is in a counter-clockwise sense. Vortex strength is T = 5.5. 

It is interesting to note that the pressure disturbances immediately downstream 
of the shock form very steep gradients. These gradients may be considered to be 
shock waves, and have been referred to as “reflected shocks” [25]. Figure 5.32 
shows the distribution of perturbation pressure as a function of distance from the 
initial shock position at locations above and below the vortex filament. It is clear 
that there are significant jumps in pressure downstream of the original shock wave, 
especially below the vortex filament. The formation of steep gradients downstream 
requires that the algorithm used in the calculation is robust. 

Figure 5.33 shows contours of the density perturbations downstream of the 
shock. The same types of features present in the T = .75 case are present here, 
but the strength of these features is now enhanced. 


144 



Distance from Core 



Distance from Core 



core at 
= 5 . 5 . 



145 



Figure 5.32: Pressure perturbations as functions of anal position at T 

Dashed curve represents pressure perturbation afong the hue r = 98, below the vor- 
tex filament . Solid curve represents pressure perturbation along the hue r = 158 

above the vortex filament. 
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Figure 5.34: Contours of pressure perturbation downstream of the shock at T — 50. 
Vortex rotation is in a clockwise sense. 

5.5 Typical Interaction-Clockwise Vortex 

In this section, results are presented for the interaction of a clock-wise rotating 
vortex with a shock. This case iB not representative of the physics of the interaction 
of vortices in the jet shear layer with shock waves, but is included for completeness. 
This type of interaction could model wake flow- shock interaction, such as the wake 
of a helicopter blade interacting with a shock on the subsequent blade. 

Figure 5.34 shows contours of pressure for the case of vortex rotating in a 
clockwise sense with a shock wave. The strength of the vortex in this case is 
T = —0.75, and the Mach number upstream of the shock is M = 1.5. 

Note that this case is directly analogous to the case presented in Section 5.3, 
except for the sign of the vortex circulation. Comparison of Figure 5.34 with Figure 
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5.4 shows that the contours of pressure perturbation are quite similar, except that 
the regions of compression (rarefaction) in Figure 5.34 are regions of rarefaction 
(compression) in Figure 5.4. The difference in the rotation sense of the vortex re- 
sults in a different response of the shock wave, which in turn results in a difference 
in the sign of the pressure disturbance downstream. As shown in Figure 5.19, a 
CCW rotating vortex causes the shock to move upstream in the region above the 
vortex filament and downstream below the vortex filament. For the clockwise ro- 
tation, the shock wave moves downstream above the vortex filament and upstream 
below the vortex filament, as shown in 5.35. Figure 5.36 clearly shows that the 
shock displacement is significantly greater upstream than downstream. 

The downstream pressure disturbances resulting from the CW and CCW vor- 
tices are not perfect images of each other about the vortex filament, the only 
difference being the sign of the disturbance. Clearly, the cusp-structure described 
in Section 5.3 connects the acoustic disturbance to a position on the shock near 
the axis of symmetry in both cases. 

There is also asymmetry of the acoustic disturbances. Figure 5.37 shows the 
pressure perturbations along radii extending from the vortex core at ±40, ±50, 
and ±60 degrees. The dashed fines represent results from the CCW vortex and 
the solid fines represent results from the clock-wise rotating vortex. For ease of 
comparison, the results have been plotted such that the pressure disturbance asso- 
ciated with the positive angle of the clockwise vortex is compared with the negative 
angle of the CCW rotating vortex. The results show good agreement between the 
pressure peak at the positive angle for the clockwise vortex and the negative angle 
for the CCW vortex. However, there is a significant difference in the maximum 
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Figure 5.35: Shock displacement as a function of space-time. 

pressure amplitude of the disturbances which correspond to the clockwise vortex 
at a negative angle and the CCW vortex at the positive angle. For these cases, 
the clockwise rotating vortex generates a larger pressure disturbance. The largest 
difference is seen to be 156 % for the sound disturbance traveling at ±60 degrees 
from the horizontal. 

5.6 Effect of Vortex Core Size 

To show the effect of vortex core size relative to the vortex ring size, a study is 
included in which the ratio of the core radius to ring radius is ^q. Figure 5.38 shows 
the contours of change in pressure at T = 50 for this case. The range of pressure 
levels in the contour plot is kept the same as in Figure 5.4 for direct comparison. 
The strength and wavelength of the acoustic waves are essentially identical to 
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Figure 5.36: Shock displacement as a function of radial distance (on the vertical 
axis) for T = 1, T = 6, and T = 10 through T = 50 in increments of 10 for a 
clock-wise rotating vortex, core at x = 30, y = 125. Results are shown for radii at 
±40, ±50 and ±60 degrees. The solid lines represent solutions for the clock- wise 
rotating vortex. The dashed lines represent solutions for the CCW rotating vortex. 
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Figure 5.37: Pressure perturbations along radii extending from the vortex cor 
x = 30, y = 125. Results are shown for radii at ±40, ±50 and ±60 degrees. The 
solid lines represent solutions for the clock-wise rotating vortex. The dashed bnes 
represent solutions for the CCW rotating vortex. 
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Figure 5.38: Contours of pressure perturbation for a case where the strength of 
the vortex is T = 0.75, the Mach number upstream of the shock is M = 1.5, and 
the ratio of the vortex core radius to the ring radius is 

those in the case where the ratio of core radius to ring radius is ^ . The primary 
difference between the results is the fact that the cusp-like structure so apparent 
in the results of Figure 5.4 is no longer visible on the contour plot. This indicates 
that the presence of this wave may be related to axisymmetry. As the vortex ring 
becomes larger relative to the core size, the resultant wave structures more closely 
resemble those generated in a two-dimensional planar interaction which does not 
show the cusp-like structure. 
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5.7 Conclusions 

In this chapter, the interaction of a ring vortex with a shock wave is presented as a 
simple model for a mechanism responsible for generating shock noise in imperfectly 
expanded supersonic jets. This model was inspired by the early works of Ribner [4] 
[6] and Moore [5] who considered two dimensional interactions of disturbances with 
shock waves in the context of linear theory, but differs from these pioneering studies 
because the interaction considered is axisymmetric, to more closely model the 
interaction of axisymmetric turbulent structures within axisymmetric jets, and the 
computational nature of the present study provides a tool for solving the nonlinear 
equations governing the fluid dynamics of this very complex interaction. 

An effort was made in the design of these calculations to model the physical 
parameters of an imperfectly expanded jet. The size of the vortex core relative 
to the vortex ring radius was small because of numerical considerations, so flow 
parameters were chosen to model the interaction of a turbulent disturbance with 
a shock wave in the first shock cell. The magnitude of the vortex strength and 
the sense of the vortex rotation were chosen to be consistent with those observed 
in experiments. The range of Mach numbers studied was within the range for 
practical nozzles. 

Observations of the evolution of perturbations in pressure, density, entropy, 
vorticity, and velocity downstream of the shock wave during an interaction are 
made. Observation of these flow quantities leads to the conclusion that sound, 
entropy, and vorticity are generated downstream of the shock wave during the 
interaction process. The sound wave is isentropic and propagates downstream at 
the sum of the convection and sound speeds. Contact surfaces are formed during 
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the interaction process. These disturbances do not support a difference in pressure, 
but do support differences in density, entropy, and vorticity. Both sound waves 
and contact surfaces have been observed in experimental studies of shock-vortex 
interaction. It is believed that these are the first calculations performed with 
enough resolution to show the presence of the contact surfaces downstream of the 
shock for these flows. Contact surfaces have been observed in experiments, which 
provides evidence for the physical nature of these disturbances, but, as discussed 
in Chapter 3, numerical error associated with moving shock waves may readdy 

manifest itself in entropy, so further analysis is necessary in order to validate the 
strength of the computed contact disturbances. 

The structure of the contact surfaces is related to the shock dynamics. The 
contact disturbances are observed to be generated by a wave which is initiated 
during the interaction and travels along the shock both towards and away from 
the center of the vortex ring. The motion of the shock corresponding to the wave 
traveling along its length generates entropy and vorticity disturbances. Pressure 
perturbations are also observed to originate at the locations of large shock motion. 

Analysis of the data provides further insight into the physics of the interaction. 

A frequency analysis suggests that the sound wave generated in these computations 
decays as a cylindrical wave. This is due to the presence of the shock wave which 
acts as a barrier to sound spreading upstream of the shock. An analysis of sound 
intensity level provides insights into the directive nature of the sound wave. High 
intensity levels are observed close to the shock wave and in downstream directions 
at angles of approximately 50 and -55 degrees from a horizontal axis through the 
vortex filament for the case where the upstream Mach number is 1.5. 
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Additional studies of the interaction of a strong vortex with a shock wave 
validated the robustness of the numerical code used in this work. The interaction 
produced extremely high gradients in the flow downstream of the shock, and the 
computation maintained stability. 

The sound generated by clockwise and counter-clockwise rotating vortices was 
compared. It was found that the structure of the alternating compression-rarefaction 
zones along the sound wave changed sign when the sense of the vortex rotation 
changed. (Regions of compression resulting from the clockwise vortex became re- 
gions of rarefaction for the counter clock- wise vortex, and vice verse.) It is also ob- 
served that the clockwise vortex generates disturbances which can be significantly 
larger in amplitude. The peak pressure perturbation resulting from the interaction 
of the clockwise rotating vortex with the shock produces pressure levels up to 156 
% higher than the counter-clockwise rotating vortex. 

A study of the effect of the ratio of core size to ring size is also performed. This 
study shows that the effects of axisymmetry are reduced when the core radius size 
is decreased relative to ring radius. 

A study of the effect of flow Mach number on the sound pressure level provided 
a significant result. It was found that the sound pressure level increased with 
Mach number, and the rate of this increase corresponds closely to that observed 
in experimental data of shock noise of supersonic jets. This result implies that the 
interaction of a vortex ring with a shock wave is a dominant physical process in 
shock noise generation in supersonic jets. 
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Chapter 6 


Conclusions 


This research shows that computational methods can predict shock-generated sound 
Although the direct computation of sound for realistic three-dimensional, time 
dependent aerodynamic flow is currently limited by computer time and memory 
constraints, if current advancements in computer hardware and software continue, 
these limitations will be removed in the near future. In the meantime, significant 
insight can be gained into the physics of sound generation in complicated flow fields 
by modeling fundamental elements of these flow fields. Understanding the sound 
generation is a first step towards the development of quieter aircraft, equipment 
and living spaces. 

This research pioneers the application of a direct computational approach to 
the study of shock noise mechanisms. Direct simulation of sound generation in 
shocked flows is challenging because of the disparity in amplitudes between the 
acoustic waves and the shocks. These challenges are met by the implementation 
of a high-order accurate Essentially Non-Oscillatory (ENO) scheme which uses 
adaptive stenciling to maintain high-order accuracy in smooth regions of the flow 
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development of sound. 

Shock motion is modeled by the interaction of a sound wave with a shock. 
Analysis of shock motion using LighthilTs equation shows that monopole, dipole, 
and quadrupole terms all have potential to contribute to the far field sound. At low 
supersonic Mach numbers, the monopole term dominates, followed by the dipole. 
The dipole term is highly directional. 

Shock motion is modeled numerically by the interaction of a sound wave with 
a shock. During the interaction, the shock wave begins to move and the sound 
pressure is amplified as the wave passes through the shock. Computations of 
sound waves interacting with shocks in a converging diverging nozzle are performed. 
The results show that the amplitude of the transmitted pressure perturbation is 
greater than the incident pressure perturbation for all Mach numbers greater than 
one. The numerical approach is validated by comparison of the computed ratio 
of transmitted to incident sound pressures with linear theory. The comparison 
is good over the range of shock Mach numbers studied (1.3 < M < 2.6). An 
energy analysis is performed to determine if acoustic energy is generated in the 
sound-shock interaction. The analysis is based on an exact representation of the 
transport of energy in an arbitrary flow field, and shows that acoustic energy is 
indeed generated during the sound-shock interaction. The source term of this 
energy is shown to be confined to a region along the shock wave in space-time, 
and is a function of the changes in entropy, momentum, and temperature from the 
mean flow state. 

Shock deformation is investigated by the simulation of a ring vortex interacting 
with a shock wave. This model has practical significance because it models the pas- 
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sage of a turbulent structure through a shock wave. Observations of the evolution 
of perturbations in pressure, density, entropy, vorticity and velocity downstream of 
the shock lead to the conclusion that acoustic waves and contact surfaces are gen- 
erated by the interaction. That these two types of fluid structures are generated by 
the interaction is validated by experimental evidence. The structure of the contact 
surfaces is related to the shock dynamics. The contact surfaces are observed to 
be generated by disturbances which are initiated during the interaction and travel 
along the shock both towards and away from the center of the vortex ring. The 
motion and deformation of the shock generates entropy and vorticity, respectively. 

Analysis of the numerical results demonstrates that the sound wave which re- 
sults from the interaction of a vortex ring with a shock wave spreads cylindrically. 
This is due primarily to the presence of the shock wave which acts as a barrier to 
sound traveling upstream. Analysis of the sound intensity level over the region of 
the computation provides insight into the directivity of the sound. High intensity 
levels are seen along the shock wave and at angles of approximately 50 and -55 de- 
grees from a horizontal axis through the vortex filament when the upstream Mach 
number is 1.5. The peak sound amplitude generated by a a clockwise rotating 
vortex was found to be as much as 156 percent higher than sound generated by 
the interaction of a counter-clockwise rotating vortex and shock wave. A signifi- 
cant result of this work is that the sound pressure level is shown to increase with 
shock strength. The relationship between the sound pressure level (SPL) and shock 
strength, defined by the parameter /? = %/ M 2 — 1, is shown to be approximately 
SPL oc /? 4 . This is consistent with experimental observations of shock noise in 
supersonic jets, and implies that the interaction of a vortex ring with a shock wave 
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is 


a dominant physical process in the physics of shock 


noise generation. 



Appendix A 


Derivation of Unsteady Shock 
Jump Relations 


In this appendix, the shock jump relations for a moving shock are derived by using 
generalized functions. The results obtained here for the continuity and momentum 
equations are also provided in [1], However, a derivation of the unsteady shock 
jump relation for the energy is not provided in [1]. 

It is important to begin by introducing the concept of a generalized function. 
Generalized functions will not be defined in a rigorous mathematical context here 
(see [1] and its references for details), but will be presented to show the connection 
between generalized and ordinary functions for clarity. 

Conventionally, a function is defined as a table of ordered pairs (z, /(z)), where 
for each z, /(z) is unique. This table may have an infinite number of ordered pairs. 
In an analogous fashion, in generalized function theory, the function /(z) is defined 
by its action on a given space of ordinary functions called test function space: 
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Figure A.l: Schematic of a discontinuous surface. 


F[<f>\ = f f{x)(j>(x)dx (A.l) 

J — OO 

where the function 4>{x) is a test function which must satisfy certain properties 
given in [1], The mapping described by Eqn. A.l is called a functional. The 
function / is now identified by the new table F[<f>],<f> £ test function space. 

It can be shown that there are an infinite number of functions <f> which satisfy the 
conditions prescribed on the space of test functions, so the table produced by Eqn. 
A.l has an uncountably infinite number of elements. The space of test functions is 
so large that the functionals on this space generated by Eqn. A.l contain not only 
ordinary functions, but additional functions. Thus, ordinary functions are a subset 
of the generalized functions. It can be shown from classical Lebesque integration 
theory that the Dirac delta function cannot be an ordinary function [1]. However, 
functions such as the Dirac function are included in the definition of generalized 
functions. Now that the concept of function space has been extended to include 
functions such as the Dirac delta, the extension of the definition of derivative is 
presented. 
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1*1 /(*, V, m) be a piecewise smooth function with one surface of discontinuity, 
-note this surface by , = 0. Such a surface is iUustrated Fig ure A.l. At this 

,Ur ^ ° f d ‘ SCOntmU,lj ’' *>“■ * j»P in the value of the function denoted by 

A/ = /(y = 0 + )-/( 9 = <T) (A2) 

Note that g - 0+ is on the side of the surface into which Vg points and that, 
following Farassat, Vg = n. 

The generalized divergence of a vector f, V-f, is related to the ordinary dive, 
gence, V . /, and the jump across the surface normal, A f, by: 


V./ = V-/+ Vg-AfS(g) 
= V • /+ n . A/%) 


(A.3) 


The Shock jump relations for a how in which a discontinuity in the flow is 

-mg may be derived by using the definition provided by E q n. A.3, because 

the differentia, forms of the laws governing the conserve., on of mass, momentum 

an energy are v^id even when discon, inuites exist within the region, as long 

as the der.vat.ves are interpreted as generahzed derivatives. The applicat.on of 

generalized derivatives to the conservation laws is most readily performed when 

the equations are in divergence form. Begm by considering the conservafon of 
mass. 


dp 

'qI + v • (pu) = o 


Applying the definition of the generahzed derivative 


(A. 4) 


, Eqn. A.3, 
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Q Q 

+ V • (/m) + A[/m] • n%) = 0 (A. 5) 

The sum of the first and third terms define the ordinary continuity equation, 

thus this sum is zero. The term j* is equal to the negative of the velocity of the 
surface, — v n . 


Thus, 


- »*)] = 0 (A. 6) 

where u n = u ■ n is the component of velocity normal to the surface g, and v n is 
the velocity of the surface. If the states upstream and downstream of a shock are 
denoted by the subscripts of 1 and 2, respectively, and the shock is not moving, 
note that Eqn. A.6 reduces to the familar Rankine-Hugoniot relation for steady 
flows in one dimension: P 2 U 2 = p\U\. 

Consider now the conservation of momentum: 


d{pui) d(puiUj) dp 
dt + dxj + d^~° 

Applying the definition of generalized derivatives, 


(A. 7) 


9jpuj) d(fntjUj) dp 

dt dxj dxi + 


A + + Ap ^ 


%) = 0 (A. 8) 


But, since §* - -v n and $*- = nj, and noting that the sum of the first three 
terms is equivalent to the well known momentum equation for continuous flows, 
Eqn. A. 8 reduces to: 


- v n ) + pm] = 0 


(A.9) 
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It is again comforting to note that when the shock is not moving, this equation 
reduces to the familar Rankine-Hugoniot relation for steady, one-dimensional flow: 
P2 U 2 + P2 = p\u\ + p : . 

Consider now the energy equation: 

Be f d{e + p) Ui n 

dt + dxi ~ ° ( A - 10 ) 

where e is the total energy: e = ph + i/m 2 - p = pH - p. The total specific 

enthalpy is represented by the symbol H. Applying the definition of the generalized 
derivative, 


Be | d(e + p)ui 
dt + 


1 Bg 


dx{ 


+ + A f( e + pK]^:)j %) = o 

Making the simplifications similar to those made for the continuity and 
mentum equations, 


(All) 


mo- 


- A[e]u n -f- A[(e + p)u n ] = 0 


(A.12) 


2 Pu2 ~ ^ Vn + ^/m 2 )u n ] = 0 


(A-13) 


A [(p^+ -pu 2 ){u n - v n ) +pv n ] = 0 (A. 14) 

For steady, one-dimensional flow, Eqn. A. 14 reduces to the familiar Rankine- 
Hugoniot relation: h 2 + \p 2 u\ = hi + \p\u\. 
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Appendix B 


Equations for the Velocity and 
Pressure of a Ring Vortex 


For the purposes of prescribing the initial condition 


for the numerical calculation, 


the vortex rmg is assumed to be iucompressible. The vortex moves relative to a 
fixed coordinate system a. a velocity equal to the sum of the mean flow velocity, 
U, and the vortex translational velocity V. A cross-section of the vortex ring is 
filus, rated in Figure B.l. The equattons for the velocity aud pressure in the flmd 

as a result o, the presence of the vortex are derived for the purposes of prescribing 
an initial condition for a numerical calculation. 

B.l Velocity 

B.l.l Outside the Core 

Lamb |1] provides the express™ for the stream function, * of a vortex rinn: 


, IYro 1 / 2 2 9 

* = ~-^~\(^-K)K(k)-l E (k)) 
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Figure B.l: Ring vortex moving at a velocity U + V with respect to a fixed coor- 
dinate system (z,r). 

where T is the circulation, r is the distance from the axis of symmetry, K and E 
are the complete elliptic integrals of the first and second kind, 

X 

K(m) = (1 — msin 2 6)~?d6 (B.2) 

and 

X 

E{m) = JJ (1 - msin 2 6)1 d6, (B.3) 

and 

£ 2 _ 4rrp 

[(z - ®o) 2 + (r + r 0 ) 2 ] 

where ro and xq are the radial and axial positions of the vortex filament, respec- 
tively. 



172 


Given the stream function, the axial and radial velocity components, u and v , 
respectively, can be found: 

„ = = + < r ° ~ ( * ~ 2 X ° )Z - ^ E(k)] (B.5) 

r dr 27T7’2 r l 

and 

v = jty = -. r (l ~ Io) [ji:(t)- (rg + (x ~f o)2+ — mi (B«) 

r dx 2 nr 2 rr 2 r i 

where n = y/(r - r 0 ) 2 + (* - zo) 2 and r 2 = y/(r + r 0 ) 2 + (* - *o) 2 - 

In the numerical implementation of these equations for the vortex velocity, the 
polynomial approximations to K(k) and E(k), which have an error bounded by 
2 x 10 -8 [2], are used. 

B.1.2 Inside the Core 

In order to avoid the mathematical singularity on the vortex filament and to better 
model the physics of a real, viscous vortex, the velocity distribution inside the 
vortex core is assumed to have a linear distribution of tangential velocity: 


u 


■ n R 

v$ sin 6 — 

r c 


(B.7) 


and 

v — vq cos 6 — (B-8) 

r c 

where v e is the tangential velocity at the core radius r c and R and 6 are the polar 
coordinates centered on the vortex filament. The tangential velocity v e is known 
by substituting (x = x 0 + r c ,r = r 0 ) into Eqn. B.6. Note that this description of 
the core assumes that the core is circular. This is an approximation since the core 
is elliptical in shape for finit values of 
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B.2 Pressure 

B.2.1 Outside the Core 

The pressure field of this vortex is determined by momentum conservation. Con- 
sider the momentum equation for an irrotational flow: 


1 q 2 du 


(B.9) 


where p is pressure, q is the magnitude of the velocity vector, and u is the velocity 
vector {u, v} T . Now, the velocity component in the axial direction is u = U +u„(x — 
((t),r) where U is the mean flow velocity, u v is the perturbation velocity induced 
by the vortex, and the coordinate £ is fixed on the vortex filament position. The 
radial component of velocity, v = v v (x - £(t),r), where v v is the velocity induced 
by the vortex. 

The flow is not steady in the fixed reference frame, because the position of 
the vortex varies as a function of time: £ = x + {U + V^f, where V is the vortex 
translational velocity. Evaluating the derivative ^ one obtains 

du _ duy{x — ^(t),y) 

9t 8u 9t 8(» - ((t)) 


d(x -£(*)) 

du v 


dt 


0 


d(x 

= - a ^ u ^ v) 


(U + V) 


(B.10) 


Similarly, 

(BU) 
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Note that for 


. i n dv flu nrViPTpfore = so that 

irrotational flow, ltierexore, ^ “5T 


— = —Vu v (U + VO 

at 


(B-12) 


If density is assumed to be constant, the momentum equation can be written: 


V[- + - — {U + V)u„] = o 


(B.13) 


where , 2 = » 2 + v 2 = (U + u„) 2 + «* Substituting in for «, v, and evaluating the 
integral of Equation B.13 along a stream line which terminates in the mean how 
state where the velocity induced by the vortex is aero, one obtains 

p - po = -^(» 2 +uS) + hoV'«v ( B -' 4 > 

where the subscript 0 denotes the mean flow state. 

B.2.2 Inside the Core 

Equation B.14 is valid for the Sow outside the vortex core, where the Sow is 
irrotational. To obtain the pressure held inside the vortex core where the Sow ,s 
rotational, consider the radial momentum equation: 


d P = P -^dr 

V 


(B-15) 


Substituting for v e and integrating from the edge of the vortex core to r 

p(r)=p(r c ) + f^ 2 -c) 


PQ v hj-J\ (B.16) 


assuming that p — po- 


The pressure p(r c ), determined from Eqn. B.14 to ensure 


continuity of pressure across 


the vortex core interface is: 

, v PO 2 

p{r c ) = Po ~ y 0 


(B.17) 
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r 




5p 5p 


Figure B.2: Profile of the pressure distribution of a counter-clockwise rotating ring 
vortex of strength T = 0.75. The figure to the right shows the pressure distribution 
over a smaller range of pressure perturbation to highlight the asymmetry in the 
pressure profile above and below the filament. 

The pressure profile of a counterclockwise vortex of strength T = 0.75 is pro- 
vided in Figure B.2. Note the asymmetry in the pressure field relative to the vortex 
filament and the region of slightly positive pressure below the filament. 

In the implementation of the vortex ring as an initial condition, the flow quan- 
tites are normalized with respect to the upstream static pressure and density, and 
the vortex core radius. 

B.3 Remarks 

Note that although the pressure and velocity are defined to be continuous at the 
rim of the vortex core, the derivatives of these quantities are not. Because the 
natural dissipation in the algorithm will automatically smooth these derivatives 
at the rim, discontinuities in derivatives at the rim are not a problem with this 
algorithm. Less robust schemes may have difficulty with this initial condition. 
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